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Abstract 

This  final  technical  report  summarizes  the  research  performed  by  the  Colorado  Center  for 
Astrodynamics  Research  (CCAR)  for  the  Naval  Research  Laboratory  (NRL)  on 
Spacebome  Differential  GPS  Applications  over  the  period  01  September  1997  through  31 
August  1999.  The  primary  focus  of  our  work  was  on  relative  state  estimation  using  GPS 
data  from  two  vehicles.  Data  from  a  spacebome  experiment  conducted  by  NASA  was 
used  for  algorithm  development  and  performance  analysis.  An  auxiliary  area  of  research 
was  the  identification  and  reduction  of  multipath  errors  in  GPS  observations.  This  effort 
is  key  to  full  utilization  of  GPS  for  high  precision  applications  such  as  attitude 
determination,  differential  GPS,  and  relative  navigation.  The  report  comprises  three 
documents  describing  the  research.  The  first  two  are  conference  papers  presented  at  ION 
GPS  99  describing  the  relative  navigation  and  multipath  research.  The  third  is  a  doctoral 
dissertation  research  proposal  written  by  Dolan  Highsmith,  summarizing  his  research  on 
this  project  and  plans  for  future  work. 
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ABSTRACT 

The  objective  of  this  research  is  to  use  flight  data  to 
develop  and  demonstrate  the  most  accurate  and  robust 
GPS-based  methods  for  relative  navigation  for  two  co¬ 
operative,  co-orbiting  vehicles  in  LEO.  A  Kalman  fil¬ 
ter  that  directly  estimates  the  relative  state  between 
two  co-orbiting  vehicles  using  single  difference  pseudo¬ 
range  observations  has  been  developed.  This  filter  is 
used  in  conjunction  with  an  existing  GPS  absolute  or¬ 
bit  determination  utility,  which  provides  an  estimate 
of  the  absolute  orbit  of  the  passive  vehicle  to  the  rela¬ 
tive  filter.  The  Clohessy-Wiltshire  equations  are  used 
for  state  and  error  propagation.  This  filter  has  been 
designed  using  flight  data  from  the  STS-80/ORFEUS- 
SPAS  flight  experiment  performed  in  December  1996. 
The  root-mean-square  (RMS)  relative  position  errors 
are  6.4  m  radial,  5.3  m  in-track,  and  1.3  m  cross-track. 
The  RMS  velocity  errors  are  0.022  m/s  radial,  0.018 
m/s  in-track,  and  0.013  m/s  cross-track.  Though  the 
relative  positions  are  no  better  than  NASA/JSC  re¬ 
sults,  this  fairly  simple  filter  achieved  an  overall  im¬ 
provement  of  more  than  70%  in  relative  velocity.  This 
is  an  excellent  beginning  upon  which  to  build  future 


adaptations,  which  may  include  relative  J2,  relative 
drag,  and  carrier  phase  processing. 


INTRODUCTION 

Accurate  relative  navigation  of  spacecraft  permits  ca¬ 
pabilities  such  as  rendezvous  and  formation  flying — 
key  aspects  of  many  current  and  future  space  mis¬ 
sions.  In  low  earth  orbit  (LEO),  such  missions  can 
use  the  Global  Positioning  System  (GPS)  signal  as 
the  primary  navigational  aid  for  relative  positioning. 
For  example,  GPS  relative  navigation  can  be  used  for 
rendezvous  operations  at  the  International  Space  Sta¬ 
tion  and  for  precise  positioning  of  distributed,  coordi¬ 
nated  remote  sensing  satellites.  Such  applications  re¬ 
quire  highly  accurate  knowledge  of  the  relative  orbits. 
However,  previous  work  in  the  application  of  GPS  to 
relative  navigation  in  space  has  not  demonstrated  on 
actual  flight  data  the  performance  required  for  precise 
applications.  Therefore,  the  objective  of  this  research 
is  to  use  flight  data  to  develop  and  demonstrate  the 
most  accurate  and  robust  GPS-based  methods  for  rel¬ 
ative  navigation  for  two  cooperative,  co-orbiting  vehi¬ 
cles  in  LEO. 

A  Kalman  filter  that  directly  estimates  the  relative 
state  between  two  co-orbiting  vehicles  using  single  dif¬ 
ference  pseudorange  observations  has  been  developed. 
This  filter  is  used  in  conjunction  with  an  existing  GPS 
absolute  orbit  determination  utility,  which  provides  an 
estimate  of  the  absolute  orbit  of  the  passive  vehicle  to 
the  relative  filter.  Use  of  single  difference  pseudor¬ 
anges  reduces  errors  inherent  in  using  GPS  by  remov¬ 
ing  common  mode  errors,  such  as  the  satellite  clock  er¬ 
ror.  Furthermore,  directly  estimating  a  relative  state 
allows  the  use  of  the  Clohessy-Wiltshire  (CW)  equa¬ 
tions  for  state  and  error  propagation.  In  fact,  in  its 
current  form  the  relative  filter  requires  no  integration. 

Specifically,  this  filter  estimates  a  relative  position 
and  velocity  and  relative  receiver  clock  bias  and  drift. 
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It  uses  the  C W  equations  with  velocity  aiding  to  prop¬ 
agate  the  position  and  velocity,  incorporating  the  ab¬ 
solute  orbit  estimate  of  the  passive  vehicle  to  compute 
the  mean  motion.  The  CW  equations  are  also  used  to 
propagate  the  covariance.  The  measurement  update 
computes  the  observed  and  predicted  single  difference 
pseudoranges,  where  the  predicted  quantity  accounts 
for  the  body-fixed  GPS  antenna  offset  from  the  center 
of  gravity  by  using  a  reported  attitude  estimate. 

Previous  GPS  relative  navigation  research  us¬ 
ing  flight  data  has  concentrated  on  determining  the 
relative  orbit  by  differencing  two  absolute  orbits. 
Schiesser  et  al.  implemented  such  a  filter  for  the 
Shuttle/ORFEUS-SPAS  relative  navigation  experi¬ 
ment  during  the  STS-80  mission  in  December  1996  [1]. 
Their  filter  uses  GPS  coarse/acquisition  (C/A)  code 
pseudoranges  to  solve  for  the  absolute  orbits  of  the 
two  vehicles,  then  differences  them  to  form  the  relative 
orbit.  Unlike  a  relative  state  filter  which  must  interpo¬ 
late  asynchronous  observations  taken  from  the  differ¬ 
ent  receivers,  the  absolute  orbit  differencing  method 
does  not  require  interpolation,  nor  does  it  necessarily 
require  observations  from  common  GPS  satellites. 

More  recently,  D ’Souza  et  al.  from  Draper  Lab 
developed  a  filter  very  similar  to  the  one  presented 
here  for  use  with  two  NASDA  projects,  Engineering 
Test  Satellite  7  (ETS-VII)  and  the  H-II  Transfer  Ve¬ 
hicle  (HTV)  [2].  Excellent  simulation  results  were  ob¬ 
tained.  Presuming  NASDA  implemented  this  filter  for 
ETS-VII,  Mokuno  et  al.  indicate  that  on-orbit  results 
were  also  good,  obtaining  a  GPS-derived  relative  range 
within  1  m  of  the  rendezvous  radar  measurement  and 
an  in-plane  velocity  error  within  2  cm/s  [3].  The  ma¬ 
jor  difference  between  the  Draper  filter  and  the  one 
developed  in  this  paper  is  the  dynamic  model  used  for 
state  propagation.  The  Draper  filter  uses  a  fourth- 
order  Runga-Kutta  integrator  with  J2,  J3,  J4,  and 
drag  to  integrate  the  absolute  states  for  both  vehi¬ 
cles,  then  creates  the  relative  state  from  their  differ¬ 
ence.  Though  undoubtedly  more  accurate  than  the 
linear,  unperturbed  CW  propagation  used  here,  the 
CW  propagation  was  found  to  be  sufficient  for  mea¬ 
surement  intervals  of  ten  seconds  or  less,  with  fast  re¬ 
coveries  after  large  data  outages. 

The  following  sections  describe  the  development  of 
the  relative  state  filter  using  the  flight  data  from  the 
STS-80  rendezvous  experiment. 


FLIGHT  EXPERIMENT 

The  STS-80  mission  launched  November  19,  1996,  and 
returned  to  Earth  December  7,  1996.  During  the  mis¬ 
sion,  the  Space  Shuttle  Columbia  released  the  free- 
flying  Orbiting  and  Retrievable  Far  and  Extreme  Ul¬ 
traviolet  Spectrograph  Shuttle  Pallet  Satellite  (OR- 
FEUS  SPAS).  Both  the  SPAS  and  the  Shuttle  had 
GPS  receivers  on  board,  which  were  used  to  record 
concurrent  observations  for  demonstration  of  relative 
navigation  using  GPS.  The  SPAS  had  a  nine-channel 
Laben  Tensor  GPS  receiver  with  four  antennas,  while 
the  Shuttle  used  a  six-channel  Trimble  Advanced  Nav¬ 
igation  System  Quadrex  GPS  receiver  with  two  anten¬ 
nas  mounted  across  the  payload  bay  [1]. 

The  data  set  used  in  this  work  comprises  the  final 
73  minutes  of  rendezvous,  up  to  the  capture  (grapple) 
of  the  SPAS  with  the  Shuttle  robotic  arm.  Several 
data  sets  were  taken  throughout  the  mission,  but  only 
the  rendezvous  portion  had  a  valid  relative  orbit  truth 
reference.  This  segment  of  the  experiment  employed  a 
laser  on  the  Shuttle  for  determining  range  and  range- 
rate  to  the  SPAS.  These  data  were  converted  into  a 
best  estimated  relative  trajectory  (BERT)  in  radial, 
in-track,  and  cross-track  (RIC)  coordinates. 

Because  of  the  frequent  thrusting  of  the  Shuttle 
during  rendezvous,  the  filter  includes  compensation  for 
the  delta-V  maneuvers.  Given  Shuttle  delta-V  vectors 
expressed  in  RIC  coordinates,  implementation  of  ve¬ 
locity  aiding  in  the  filter  is  straightforward.  The  delta- 
V  measurements  were  produced  by  the  Shuttle  IMU 
accelerometer,  which  reported  accumulated  sensed  ve¬ 
locity  [1]. 

Since  the  relative  navigation  filter  relies  on  sin¬ 
gle  difference  pseudoranges,  only  common-view  GPS 
satellite  observations  may  be  used.  Fortunately,  the 
rendezvous  phase  of  the  flight  was  performed  in  such 
as  way  as  to  maximize  the  potential  for  common  view 
satellites.  The  Shuttle  payload  bay  remained  pointed 
at  the  SPAS  so  the  crew  could  maintain  visual  contact; 
consequently,  the  SPAS  was  commanded  to  pitch  at  a 
constant  rate  so  that  its  GPS  antennas  pointed  in  ap¬ 
proximately  the  same  direction  as  those  on  the  Shuttle 
throughout  the  rendezvous  [1]. 


OBSERVATIONS 

Both  GPS  receivers  reported  C/A  code  phase,  LI  car¬ 
rier  phase,  and  LI  Doppler.  The  filter  described  in 
this  paper  uses  the  C/A  code  pseudoranges  from  each 
common- view  GPS  satellite  to  create  the  single  differ- 
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ence  observation.  The  receivers  were  not  synchronized, 
so  the  measurement  time  tags  are  different.  Therefore, 
the  pseudoranges  from  the  active  vehicle  (Shuttle)  re¬ 
ceiver  are  interpolated  to  the  times  of  the  passive  ve¬ 
hicle  (SPAS)  pseudoranges  using  polint.c,  a  polyno¬ 
mial  interpolation  routine  from  Numerical  Recipes  [4]. 
Polint.c  uses  Neville’s  algorithm  to  construct  the  in¬ 
terpolating  polynomial.  As  implemented  in  the  filter, 
the  desired  point  is  computed  between  the  second  and 
third  values  of  a  four-point  interpolation. 

Given  a  noise  level  on  each  C/A  code  pseudor¬ 
ange  of  approximately  3  meters,  interpolated  points 
between  four  closely-timed  observations  will  also  be 
noisy.  This  effect  contributes  to  large  measurement 
residuals  obtained  using  these  data.  A  potential 
remedy  for  the  large  noise,  short  of  using  carrier 
phase  pseudoranges,  is  carrier  interpolation  or  carrier 
smoothing  of  the  C/A  code  pseudoranges.  Preliminary 
results  using  these  methods,  however,  do  not  show  suf¬ 
ficient  improvement  to  warrant  inclusion  in  this  paper. 


RELATIVE  STATE  KALMAN  FILTER 


Overview 

The  relative  state  sr  consists  of  a  relative  position  rr 
and  velocity  vr,  as  well  as  relative  receiver  clock  bias 
br  and  frequency  fr,  such  that 

8r=[r?  ,r  br  frf.  (1) 

The  relative  state  is  referenced  to  the  passive  vehi¬ 
cle  (SPAS).  The  filter  position  and  velocity  coordinate 
frame  is  centered  on  the  passive  vehicle  and  expressed 
in  radial,  in-track,  and  cross-track  directions.  Like¬ 
wise,  the  covariance  and  process  noise  values  are  spec¬ 
ified  for  the  RIC  directions.  This  formulation  permits 
more  direct  use  of  the  CW  equations  for  state  and  co- 
variance  propagation. 

Input  data  to  the  Relative  State  Kalman  Filter 
(RSKF)  include  asynchronous  GPS  code  pseudoranges 
for  both  vehicles,  absolute  orbit  estimates  for  the  pas¬ 
sive  vehicle,  and  GPS  satellite  position  vectors.  Since 
the  RSKF  is  designed  to  work  in  tandem  with  an  ex¬ 
ternal  GPS  absolute  orbit  determination  (AOD)  util¬ 
ity,  output  from  it  is  used  as  input  to  the  RSKF .  The 
external  AOD  program  generates  the  absolute  orbit 
estimates  of  the  passive  vehicle  using  GPS  code  pseu¬ 
doranges  and  broadcast  GPS  satellite  orbits.  The  orbit 
estimates,  along  with  the  GPS  observations  and  satel¬ 
lite  positions  used  to  compute  them,  are  output  from 
the  AOD  program  for  use  as  input  to  the  RSKF.  As 


mentioned  previously,  within  the  RSKF  the  active  ve¬ 
hicle  pseudoranges  are  interpolated  to  the  times  of  the 
passive  vehicle  pseudoranges  and  positions. 

Inertial-to-body  attitude  quaternions  and  delta-V 
vectors  are  also  input  to  the  RSKF.  The  attitude 
quaternions  of  the  Shuttle  are  needed  to  account  for 
the  GPS  antenna  CG  offsets,  which  are  specified  in  a 
body-fixed  frame.  Since  the  SPAS  attitude  mimicked 
that  of  the  Shuttle,  a  total  relative  offset  is  computed 
in  the  Shuttle  body  frame,  which  is  then  rotated  to 
RIC  using  the  Shuttle  attitude.  The  delta-V  vectors 
are  used  in  the  time  update  propagation  to  account 
for  the  thrusting  of  the  Shuttle  during  rendezvous. 


Time  Update 


As  stated  previously,  the  CW  equations  are  used  in 
this  filter  for  relative  state  and  covariance  propagation. 
From  Prussing  and  Conway  [5],  the  CW  equations  in 
state  transition  matrix  form  are 


*r,v  (A t)  = 


'  4 -3c 

0 

0 

TL 

2(1— c) 

fl 

0  ' 

6(s  —  nAt) 

1 

0 

2(1  -c) 
n 

4s— 3nAt 
n 

0 

0 

0 

c 

0 

0 

8_ 

n 

3  ns 

0 

0 

c 

2s 

0 

— 6n(l  —  c) 

0 

0 

-2s 

4c -3 

0 

0 

0 

—ns 

0 

0 

c 

where  n  signifies  the  mean  motion 

of  the  passive 

hide,  At  is  the  time  interval  for  propagation,  and 
s  =  sin(nAi)  and  c  —  cos(nAt).  In  the  filter,  the 
mean  motion  in  the  equation  above  is  actually  com¬ 
puted  as  the  angular  rotation  rate  (6)  at  the  current 
time.  In  other  words,  the  mean  motion  is  chosen  at 
each  epoch  to  be  the  current  along-track  velocity  of  the 
passive  vehicle  divided  by  its  radius.  This  provides  a 
more  accurate  propagation  for  elliptical  orbits  as  long 
as  the  measurement  interval  is  not  large. 


The  relative  clock  frequency  is  assumed  constant 
over  the  time  interval,  so  the  clock  terms  are  propa¬ 
gated  with  [6] 


$(,,/  (At)  = 


1  At 
0  1 


(3) 


Thus,  the  full  relative  state  transition  matrix  has  the 
form 


$  (At)  = 


$r,v  (At)  06x2  1 

02x6  $b,f  (At)  J 


(4) 


Incorporated  into  the  state  propagation  is  velocity 
aiding,  accomplished  by  using  the  Av  vectors  of  the 
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active  vehicle.  The  velocity  vectors  are  added  directly 
into  the  linear  state  propagation  at  the  given  time  as 
impulsive  maneuvers,  with  potentially  more  than  one 
maneuver  between  measurement  epochs.  The  follow¬ 
ing  example  illustrates  the  equations  for  the  occur¬ 
rence  of  one  Av  during  the  interval.  As  such,  the  state 
must  be  propagated  from  tk  to  tf~+i  with  a  maneuver 
at  tAv  such  that  tk  <  tAv  <  4+i*  In  this  case,  the 
state  vector  is  propagated  as  follows: 


SA»  =  $(*a®,4)  si 

(5) 

£.  +  [<&!  Av^oLif 

(6) 

Sfc+1  =  *(4+i,*At>)  s£„ 

(7) 

where  the  “  and  +  indicate  state  vector  values  before 
or  after  the  measurement  update  or  a  maneuver.  The 
distinction  is  clear  from  the  context. 

Not  portrayed  in  the  above  equations  is  the  fact 
that  the  CW  state  propagation  is  performed  in  cylin¬ 
drical  coordinates,  with  the  along-track  direction  be¬ 
ing  along  the  arc  of  the  orbit.  Since  this  filter  oper¬ 
ates  in  the  rectangular  RIC  coordinates,  a  rotation  is 
performed  to  cylindrical  RIC  coordinates  immediately 
before  each  state  propagation  step  (e.g.,  Eq.  5),  then 
back  to  rectangular  immediately  after.  The  equations 
for  such  a  conversion  are  given  in  [8].  However,  for  this 
filter,  those  equations  are  re-derived  for  the  case  that 
r  is  not  zero,  providing  a  slight  increase  in  accuracy. 

Ignoring  any  potential  maneuvers,  the  covariance  is 
propagated  over  the  entire  measurement  interval  using 

Pfc+1  =* (4+i, 4)  Pjt*(4+i,4)T  +  Qfc.  (8) 

The  process  noise  matrix  Q  is  derived  from  Gelb,  equa¬ 
tion  3.6-13  [7],  resulting  in  the  equation 

tfc+i 

Qjfc  =  J  $  (4+1 ,  t)  Q  (r)  $T  (4+1 ,  r)  dr.  (9) 

tk 

where  Q (t)  is  defined  to  be  the  spectral  density  matrix. 
Eq.  9  is  solved  symbolically  using  Eq.  4  and  a  diagonal 
representation  of  the  spectral  density  matrix  Q  (t). 


Measurement  Update 


The  predicted  single  difference  pseudorange  Psd  is 
computed  from  the  geometric  relative  range,  relative 
clock  bias  estimate,  and  antenna  CG  offsets.  Because 
the  pseudoranges  in  the  above  equation  are  measured 
from  the  GPS  antennas,  psd  is  not  referenced  to  the 
centers  of  gravity  of  the  two  vehicles.  Therefore,  the 
prediction  of  the  single  difference  pseudorange  must 
account  for  the  antenna  CG  offset  vector  tant •  In 
the  equations  below,  tant  can  be  viewed  as  the  sum 
of  the  antenna  offset  vectors  from  the  two  vehicles. 


The  geometric  relative  range  is  found  by  differenc¬ 
ing  the  ranges  to  the  active  and  passive  vehicles  us¬ 
ing  the  GPS  satellite  position  Rgfs,  passive  vehicle 
absolute  position  estimate  r p,  and  predicted  relative 
position  estimate  trel  from  the  time  update.  The 
computed  ranges  to  the  active  and  passive  vehicle  re¬ 
ceivers  are 

fp  =  |Rgps  -  rp|  (11) 
?A  —  |Rgps  “  (rp  +  *rel  +  r^Arr)!  •  (12) 

The  predicted  single  difference  is  then  computed  from 

Psd  =  rA-rP  +  bREL-  (13) 


The  single  difference  pseudorange  residual  is  formed 
using 

z  =  psd- Psd-  (14) 

The  partial  differential  of  Eq.  13  with  respect  to  the 
relative  state  (Eq.  1)  produces  the  measurement  gra¬ 
dient  matrix  H  with  the  first  three  elements  being 


#n  = 
#12  = 


XgPS  “  (xp  +  Xrel  +  Xant) 


ta 


Ygps  -  ( yp  +  Vrel  +  Vant) 


ta 


#13  =  “ 


ZgPS  -  ( Zp  +  ZREL  +  zAnt) 


Ta 


so  that  the  full  matrix  is 


(15) 

(16) 
(17) 


H  =  [  Hu  #12  #13  o  0  0  1  0  ]  .  (18) 


The  Kalman  gain  is  computed  from  the  equation 


K  =  Pfc"+1HT  (HP*+1HT  +  R)"1  (19) 


The  RSKF  performs  a  measurement  update  for  each 
GPS  observation  separately,  so  the  equations  in  this 
section  refer  to  quantities  for  one  common  view  GPS 
satellite.  Given  pseudorange  observations  to  the  same 
GPS  satellite  from  the  active  receiver  pa  and  passive 
receiver  pp ,  the  single  difference  observation  is  simply 

Psd  =  Pa~Pp-  (10) 


and  the  covariance  is  updated  with 

P++1  =  (I  -  KH)  P^+1  (I  -  KH)r  +  KRKt.  (20) 

Since  the  measurement  update  is  performed  for  each 
GPS  satellite  separately,  the  measurement  noise  ma¬ 
trix  R  is  1  x  1,  so  the  inversion  in  Eq.  19  is  a  simple 
scalar  inversion. 
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Table  1:  Relative  orbit  position  error  statistics  with 
respect  to  laser  BERT  for  the  Relative  State  Kalman 
Filter  and  the  NAS  A/ JSC  filter  [1]. 


Source 

m 

Rad 

In 

Cross 

Mean 

-2.8 

-4.2 

0.2 

RSKF 

STD 

5.8 

3.2 

1.3 

RMS 

6.4 

5.3 

1.3 

RSS 

8.4 

Mean 

-1.5 

-2.6 

-0.2 

Ten 

STD 

7.5 

2.9 

2.9 

RMS 

7.6 

3.9 

2.9 

RSS 

9.0 

Table  2:  Relative  orbit  velocity  error  statistics  with 
respect  to  laser  BERT  for  the  Relative  State  Kalman 
Filter  and  the  NAS  A/ JSC  filter  [1]. 


Source 

cm/s 

Rad 

In 

Cross 

Mean 

-1.0 

-0.4 

0.1 

RSKF 

STD 

2.0 

1.8 

1.3 

RMS 

2.2 

1.8 

1.3 

RSS 

3.1 

Mean 

-2.0 

0.0 

-1.2 

TCP 

STD 

8.7 

6.1 

4.0 

JuVj 

RMS 

8.9 

6.1 

4.2 

RSS 

11.6 

Finally,  the  relative  state  is  updated  by 

si+i  =  sk+i+Kz-  (2i) 

The  measurement  update  is  repeated  for  each  common 
view  GPS  satellite.  After  exhausting  all  the  observa¬ 
tions  for  this  epoch,  the  next  measurements  are  chosen 
and  the  entire  process  is  repeated. 


RESULTS  &  DISCUSSION 

Given  noisy,  asynchronous  observations  from  the  drag¬ 
laden,  maneuvering  Shuttle  and  the  relatively  well- 
behaved  SPAS,  the  RSKF  obtains  excellent  results 
considering  its  simplicity.  In  fact,  the  results  are 
comparable  or  better  than  results  obtained  in  a 
NASA/JSC  study  detailed  in  [1].  Tables  1  and  2  show 
the  position  and  velocity  error  statistics  for  the  same 
73  minute  period  using  the  same  truth  reference  (the 
laser  BERT)  for  the  RSKF  developed  in  this  paper 
and  the  filter  developed  by  JSC.  These  results  are  ob¬ 
tained  with  the  following  initial  variances:  4000,  1000, 


Figure  1:  In-plane  trajectories  for  filter  and  laser 
BERT. 


and  100  m2  for  RIC  relative  position;  0.004,  0.004, 
and  0.001  m2/s2  for  RIC  relative  velocity;  2000  m2 
for  relative  clock  bias;  and  25  m2/s2  for  relative  clock 
frequency.  The  initial  spectral  densities  are  as  follows: 
le-5,  4e-6,  8e-6  m2/s  for  RIC  relative  position;  3e-7, 
4e-8,  8e-8  m2/s3  for  RIC  relative  velocity;  0.01  m2/s 
for  relative  clock  bias;  and  0.008  m2/s3  for  relative 
clock  frequency.  The  noise  on  the  measurement  is  es¬ 
timated  to  be  approximately  3  m,  so  the  scalar  mea¬ 
surement  noise  value  R  is  set  to  9  m2.  It  becomes 
apparent  in  subsequent  plots  that  additional  tuning 
must  be  performed  so  that  the  filter  3 a  RMS  values 
more  accurately  reflect  the  behavior  of  the  system. 

To  illustrate  the  nature  of  the  rendezvous,  Figs.  1, 
2,  and  3  show  the  relative  position  and  velocity  states. 
The  in-plane  relative  trajectories  for  the  filter  and  laser 
BERT  are  shown  in  Fig.  1.  Note  that  motion  is  to  the 
right  and  the  Earth  is  down.  It  can  be  seen  from  this 
plot  that  the  rendezvous  data  set  begins  at  a  range 
of  approximately  1.8  km  and  ends  at  SPAS  capture. 
Figs.  2  and  3  show  all  the  relative  position  and  ve¬ 
locity  states  as  determined  by  the  RSKF  and  the  laser 
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Figure  2:  Radial,  along-track,  and  cross-track  relative 
positions  for  the  filter  solutions  and  laser  BERT. 


BERT.  These  plots  reveal  the  large  dynamic  motion  of 
the  Shuttle  during  the  first  fifteen  minutes  as  it  nar¬ 
rows  the  relative  range  from  1.8  km  to  200  m.  After 
the  first  fifteen  minutes,  the  relative  motion  becomes 
fairly  benign. 

Another  indicator  of  the  large  motions  is  the  plot 
of  the  delta-V  magnitudes  in  Fig.  4.  For  the  purpose 
of  velocity  aiding  in  the  RSKF,  a  range  of  minimum 
thresholds  was  evaluated.  A  threshold  of  1.5  cm/s 
was  determined  to  provide  the  best  accuracy.  The  fig¬ 
ure  shows  a  line  at  1.5  cm/s,  with  the  points  above 
it  marked  by  squares  to  indicate  those  used  in  state 
propagation. 

The  relative  position  errors  with  respect  to  the  laser 
BERT  in  the  radial,  along-track,  and  cross-track  direc¬ 
tions  are  shown  in  Fig.  5.  The  3<r  bounds  determined 
by  the  RSKF  variances  are  shown  as  the  dashed  lines. 
Fig.  6  shows  the  relative  velocity  errors  with  respect 
to  the  laser  BERT.  Again,  the  3 a  bounds  determined 
by  the  filter  variances  are  shown  as  the  dashed  lines. 
Compare  the  behavior  of  these  error  plots  with  the  mo¬ 


Figure  3:  Radial,  along-track,  and  cross-track  relative 
velocities  for  the  filter  solution  and  laser  BERT. 

tion  shown  in  Figs.  2  and  3.  Clearly,  the  RSKF  has 
difficulty  obtaining  accurate  results  during  the  highly 
dynamic  motion  of  the  first  fifteen  minutes,  even  with 
velocity  aiding.  After  that  period,  however,  the  filter 
estimates  become  less  erratic. 

Notice  the  biases  and  systematic  variations  in  the 
relative  position  error  for  the  radial  and  along-track 
directions.  These  are  reflected  in  the  mean  and  RMS 
statistics  in  Table  1.  We  are  investigating  possi¬ 
ble  sources  of  error,  including  inaccurate  GPS  an¬ 
tenna  CG  offsets,  time  tagging  errors,  or  simply  un¬ 
known/unmodeled  biases.  The  fact  remains,  however, 
that  these  biases  can  be  removed  by  adjusting  the 
body-fixed  antenna  CG  offsets  by  the  bias  amount. 
Doing  so  lowers  the  position  RMS  errors  to  the  same 
level  as  the  standard  deviations,  thereby  reducing  the 
total  RSS  position  error  from  8.4  m  to  6.9  m.  There  is 
only  a  minimal  improvement  in  the  velocity  statistics. 

To  ensure  that  the  RSKF  CW  propagation  is  not 
introducing  a  bias  to  the  results,  a  hardware-in- the- 
loop  simulation  data  set  used  in  [9]  was  used  to  verify 
its  accuracy.  A  5  km  along-track  separation  was  simu- 
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Figure  4:  Shuttle  IMU-sensed  AV  magnitudes. 
Threshold  of  1.5  cm/s  is  shown,  with  squares  signi¬ 
fying  the  points  used  in  state  propagation. 


lated,  with  the  receiver  reporting  P-code  pseudoranges 
without  selective  availability  or  multipath.  For  this 
low-noise  (approximately  3  cm)  data  set,  the  RSKF 
obtained  an  unbiased  relative  position  RSS  of  0.18  m 
and  a  relative  velocity  RSS  of  0.5  cm/s.  Clearly,  the 
source  of  the  bias  in  the  Shuttle/ SPAS  data  set  re¬ 
quires  further  investigation. 

The  most  dramatic  result  of  this  research  is  the  im¬ 
provement  in  relative  velocity  estimates  by  the  RSKF 
versus  the  JSC  results.  The  RSKF  total  RSS  velocity 
error  of  3.1  cm/s  is  a  73%  improvement  over  the  11.6 
cm/s  value  from  JSC  (see  Table  2).  Both  filters  include 
velocity  aiding  using  the  Shuttle  IMU  accelerometer 
data.  In  fact,  both  filters  use  virtually  the  same  data, 
so  the  improvement  in  velocity  accuracy  is  most  likely 
an  artifact  of  estimating  a  relative  state  as  opposed  to 
two  absolute  states. 

The  error  in  relative  semimajor  axis  (SMA)  is 
shown  in  Fig.  7.  This  error  is  computed  as  the  differ¬ 
ence  between  the  relative  SMA  computed  with  RSKF 


Figure  5:  Relative  position  error  in  radial,  along-track, 
and  cross-track  directions  with  respect  to  laser.  Dot¬ 
ted  lines  show  filter  3cr  RMS  limits 

states  and  the  relative  SMA  computed  with  the  laser 
BERT.  Note  that  for  lack  of  a  valid  absolute  orbit 
truth  for  the  passive  vehicle,  the  relative  SMA  from 
the  BERT  is  computed  using  the  passive  vehicle  ab¬ 
solute  state  from  the  external  AOD  filter.  Thus,  the 
relative  SMA  comparison  does  not  include  errors  due 
to  inaccuracies  in  the  absolute  orbit  estimate  of  the 
SPAS. 


CONCLUSIONS 

In  a  relatively  simple  form,  the  RSKF  has  demon¬ 
strated  on  flight  data  the  ability  to  predict  relative 
position  and  velocity  to  accuracy  levels  comparable  to 
or  better  than  the  absolute  state  filter  developed  by 
NASA  /JSC.  Further  improvements  in  positioning  ac¬ 
curacy  may  be  obtained  by  identifying  the  source  of 
the  relative  position  estimate  biases,  which  appear  in 
both  the  RSKF  and  the  JSC  filter  results  to  differ¬ 
ing  degrees.  There  is  also  the  potential  for  additional 
overall  improvements  in  accuracy  by  implementing  a 
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Time  (min) 


Figure  6:  Relative  velocity  error  in  radial,  along-track, 
and  cross-track  directions  with  respect  to  laser.  Dot¬ 
ted  lines  show  filter  3 <j  RMS  limits 


relative  state  propagator  that  accounts  for  the  effects 
of  relative  gravity  and  drag  perturbations.  The  most 
dramatic  improvement,  however,  is  gained  by  using 
more  precise  observations  and  techniques,  namely  car¬ 
rier  phase  processing  coupled  with  synchronization  of 
the  vehicle  clocks.  This  is  the  direction  of  future  re¬ 
search  on  the  RSKF. 
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ABSTRACT 

A  process  was  developed  to  use  the  GPS  signal  to  noise 
ratio  (SNR)  to  identify  an  effective  reflector  for  a  given 
antenna.  This  is  termed  an  effective  reflector  because  the 
reflector  identified  can  account  for  a  number  of  effects 
found  in  the  signal,  including  such  things  as  multiple 
reflectors,  changes  in  reflectivity,  and  changes  in  signal 
phase  upon  reflection.  A  physical  location  is  used  simply 
as  a  means  to  represent  the  spatial  correlation  of  the 
multipath.  Once  a  effective  reflector  location  has  been 
found,  a  carrier  phase  correction  profile  is  computed  for 
that  reflector.  The  measured  carrier  phase  is  then 
corrected  by  that  amount  and  a  new  attitude  solution  is 
computed. 

Using  this  technique,  results  from  ground  tests  show  that 
the  RMS  differential  phase  residuals  improved  from  1 1.37 
mm  to  8.56  mm,  corresponding  to  a  25%  improvement. 
These  results  were  found  using  only  a  batch  estimator  to 
identify  an  effective  reflector.  For  the  CRISTA-SPAS 
flight  data  the  lack  of  knowledge  of  the  antenna  gain 
pattern  limited  the  ability  of  the  method  to  determine  a 
reliable  reflector  solution.  Some  residual  improvement 


was  found  for  the  data.  One  such  pass  had  a  residual 
improvement  from  6.4  mm  to  5.5  mm.  However,  for 
much  of  the  data,  an  effective  reflector  that  improved  the 
phase  residuals  could  not  be  identified. 

INTRODUCTION 

High  precision  GPS  applications  such  as  attitude 
determination,  static  and  kinematic  survey  utilize  carrier 
phase  data  from  two  or  more  GPS  antennas  and  receivers. 
While  this  technique  is  very  effective  in  getting  results  at 
the  several  centimeter  level,  the  factor  limiting  accuracy 
for  short  antenna  separations  is  often  multipath. 

Much  work  has  been  done  in  the  area  of  GPS  multipath 
identification  and  correction.  One  approach  is  to  correct 
for  multipath  within  the  receiver,  as  opposed  to  a  post¬ 
processing  technique.  A  few  of  the  proposed  internal 
receiver  techniques  include  narrow  correlators  [Van 
Dierendonck,  et  al.,  1992],  edge  correlators,  and  the 
multipath  estimating  delay  lock  loop  (MEDLL) 
[Townsend,  et  al.,  1995].  These  methods  can  potentially 
reduce  the  amount  of  multipath  in  the  observable  at  the 
expense  of  increased  receiver  complexity. 

Another  approach  is  to  model  the  multipath  for  a  specific 
antenna  and  environment  using  the  geometrical  theory  of 
diffraction  (GTD)  [Hajj,  1990  and  Lippencott,  et  al., 
1993].  This  approach  requires  a  knowledge  of  the 
antenna  gain  pattern  and  the  geometry  and  reflectivity  of 
objects  surrounding  the  antenna.  Based  on  this 
information,  a  ray  tracing  algorithm  determines  a  model 
of  the  resulting  phase  and  amplitude  of  the  multipath 
error.  Corrections  on  the  carrier  have  been  confirmed  by 
experimental  results  [Gomez,  et  al.,  1995  and  Irish,  et  al., 
1998]. 

A  method  for  post-processing  of  the  GPS  carrier  phase 
signal  using  a  procedure  based  upon  the  spatial 
correlation  characteristics  of  multipath  was  first  proposed 
by  Cohen  and  Parkinson  [1991].  Experimental  data 
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collected  on  a  fixed  ground  structure  over  the  course  of 
several  days  is  used  to  create  a  map  of  multipath  as  a 
function  of  the  direction  of  the  incoming  direct  GPS 
signal.  Based  on  the  known  attitude  and  baseline 
locations  of  the  ground  platform,  carrier  phase  residuals 
are  computed,  and  prefiltered  to  remove  receiver  noise. 
Then,  an  eighth  order  spherical  harmonic  model  is  fit  to 
the  residuals  as  a  function  of  azimuth  and  elevation  of  the 
incoming  signals.  Reductions  of  the  residuals  from  5.2  to 
3.2  millimeters  were  shown  by  Cohen  and  Parkinson 
[1991]  using  this  approach. 

Another  method  for  the  post-processing  the  GPS  carrier 
phase  was  implemented  by  Comp  [1996].  Instead  of 
modeling  multipath  in  the  residual  phase,  this  approach 
utilizes  the  signal  to  noise  ratio  (SNR)  in  order  to  mitigate 
multipath.  One  advantage  of  using  the  SNR  over  the 
residual  differential  phase,  as  implemented  by  Cohen  and 
Parkinson  [1991],  is  that  the  SNR  is  much  less  sensitive 
to  attitude  errors  than  the  residual  differential  phase, 
which  is  directly  dependent  on  any  errors  in  the  attitude. 
The  previous  SNR  method  analyzes  the  data  for 
individual  satellite  passes  as  a  function  of  time, 
identifying  individual  frequency  components  in  the  signal 
due  to  multipath.  Although  this  approach  showed 
significant  reduction  in  the  phase  multipath  errors,  it  is 
not  entirely  effective  for  real-time  applications  because  of 
the  long  convergence  time  required. 

Ray,  et  al.  [1998]  developed  a  method  to  reduce  multipath 
by  using  a  system  of  multiple,  closely  spaced  antennas  in 
a  static  multipath  environment.  With  the  location  of  the 
antennas  precisely  known,  it  is  possible  to  isolate  the 
direct  signal  phase  and  discard  the  portion  of  the  signal 
corrupted  by  multipath.  A  drawback  of  this  system  is  that 
it  relies  on  using  the  phase  measurement,  which  is  highly 
sensitive  to  errors  in  antenna  location  and  attitude.  Also, 
it  relies  on  using  a  number  of  closely  spaced  antennas. 
Although  it  may  be  possible  to  do  this  for  some  ground 
stations,  it  would  be  too  costly  to  use  in  a  spacecraft 
environment  for  attitude  applications. 

The  purpose  of  this  research  is  to  avoid  the  approaches 
that  increase  the  cost  or  complexity  of  an  experiment  by 
altering  the  receiver  architecture  or  including  additional 
antennas.  A  method  will  be  implemented  in  which  post¬ 
processing  data  from  any  type  of  receiver-antenna  setup 
can  be  used.  Additionally,  the  GPS  SNR  measurement  is 
desirable  for  use  in  a  post-processing  technique,  as  was 
used  in  the  research  by  Comp  [1996],  because  the  SNR  is 
less  sensitive  than  the  phase  residuals  to  errors  in  the 
attitude  and  baseline  solutions. 


This  paper  discusses  a  method  by  which  to  reduce  the 
multipath  in  the  carrier  phase  by  utilizing  the  spatially 
correlated  characteristics  of  multipath,  as  viewed  in  both 
the  GPS  carrier  phase  and  SNR.  The  spatial  dependence 
of  multipath  is  characterized  by  an  effective  reflector 
location  and  corresponding  multipath  signal  amplitude. 
Therefore,  when  an  effective  reflector  is  determined  for  a 
set  of  SNR  data  at  some  time  from  a  certain  portion  of  the 
sky,  the  same  reflector  can  be  used  for  data  at  a  later  point 
in  time  if  the  originating  signal  passes  through  the  same 
region  of  the  sky.  The  location  and  amplitude  of  the 
effective  reflector  are  in  turn  used  to  generate  corrections 
in  the  carrier  phase  that  are  caused  by  the  effective 
reflector. 

In  the  next  section  a  discussion  of  carrier  phase  multipath 
is  given,  followed  by  the  methodology  behind  the 
reflector  identification  algorithm.  Finally,  results  are 
presented  for  static  ground  data  and  for  flight  data  from 
the  CRISTA-SPAS  spacecraft. 

BACKGROUND 

Multipath  occurs  when  an  electromagnetic  signal  arrives 
at  an  antenna,  not  along  a  direct  path,  but  via  one  of  the 
multiple  paths  the  signal  takes  as  it  reflects  off  an  object 
near  an  antenna.  Multipath  may  either  be  diffuse  or 
specular  in  nature.  In  contrast  to  the  correlated  nature  of 
specular  multipath,  diffuse  multipath  takes  on  an 
unbiased,  random  appearance,  and  is  removed  more  easily 
in  filtering  than  specular  multipath.  Therefore,  only 
specular  multipath  will  be  considered  in  this  research. 

Path  Delay 

When  a  signal  enters  the  antenna  via  an  indirect  route,  the 
extra  distance  the  signal  travels  as  compared  to  a  direct 
line  of  sight  to  the  satellite  is  defined  as  path  delay.  For 
the  purposes  of  this  study,  we  model  only  a  perfectly 
conducting  planar  reflector,  for  which  the  angle  of 
reflection  is  equal  to  the  angle  of  incidence. 

The  path  delay  due  to  a  planar  reflector  is  illustrated  in 
Figure  1.  In  this  figure,  h  is  the  normal  vector  of  the 
plane,  d  is  the  normal  distance  from  the  plane  reflector  to 
the  antenna,  and  e  is  the  unit  line  of  sight  vector  from  the 
user  to  the  GPS  satellite.  The  physical  path  delay  is 
shown  in  the  figure  as  the  bold  part  of  the  reflected  signal. 
An  expression  for  the  path  delay,  V|/,  is  given  in  equation 
(1).  Note  in  equation  (1)  that  the  path  delay  is  only 
modeled  as  a  function  of  the  distance  of  the  reflector  from 
the  antenna  (< d ),  the  orientation  of  the  reflector  («),  and 
the  direction  of  the  line  of  sight  vector  ( e ),  thus  making  it 
a  purely  spatially  correlated  phenomenon. 
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Figure  1 :  Path  Delay  for  a  Reflected  Signal 

\|f  =2d(eh)  (1) 

Multipath  in  the  GPS  SNR  and  Differential  Phase 

The  path  delay  due  to  a  reflector  is  identifiable  in  three 
measurements  from  a  GPS  receiver  -  the  pseudorange, 
the  carrier  phase  and  the  SNR.  Here  we  focus  on  the 
latter  two  measurements  depicted  in  the  phasor  diagram 
shown  in  Figure  2.  The  amplitude  values,  denoted  by  A, 
correspond  to  SNR  measurements  and  the  angular  values, 
denoted  by  <j),  correspond  to  phase  measurements.  In  this 
figure  Ac  is  the  composite  SNR  measurement  and  §c  is  the 
composite  carrier  phase.  These  are  the  quantities  that  the 
receiver  measures.  In  addition  to  the  composite  signal 
phasor,  two  others  are  represented  in  this  figure  -  the 
direct  signal  phasor  and  the  multipath  phasor. 


The  direct  signal  SNR  amplitude  is  AaA0 ,  where  Aa  is  the 
normalized  antenna  gain  and  A0  is  the  constant  part  of  the 
direct  signal  amplitude.  The  multipath  amplitude  is  Am 
and  the  path  delay,  \| /,  is  the  phase  of  the  multipath 
relative  to  the  direct  signal.  Using  the  law  of  cosines,  the 
equation  for  the  composite  SNR  in  terms  of  the  direct  and 
multipath  signals  is  shown  in  equation  (2). 

=  (A0Aa  f+A2m  +  2 A0  Aa  Am  cos X?  (2) 


The  composite  carrier  phase  measurement,  <j>c,  is  also  a 
function  of  the  multipath  amplitude  and  phase.  The  phase 
due  to  multipath  is  the  residual  phase,  &|> .  The  residual 

phase  is  represented  using  equation  (3). 


tan  (&)))= 


K  sin  V 

AoAa  +  Am  COSV 


0) 


In  order  to  gain  a  better  understanding  of  the  dynamics  of 
this  problem  note  that  as  a  GPS  satellite  passes  overhead, 
the  relationship  between  the  line  of  sight  vector  and  the 
vector  normal  to  the  plane  will  change,  causing  the  path 
delay  to  change.  As  the  path  delay  increases  or  decreases, 
the  multipath  phasor  rotates  about  the  end  point  of  the 
direct  phasor.  This  in  turn  causes  the  composite  SNR  and 
carrier  phase  magnitudes  to  oscillate  sinusoidally  as  the 
path  delay  either  increases  or  decreases,  thus  creating  the 
correlated  error  found  in  the  two  signals. 

Reflector  Identification  Algorithm 

A  batch  estimation  filter  is  implemented  in  order  to 
identify  an  effective  reflector  location.  The  batch 
algorithm  utilizes  the  dependence  on  the  reflector  location 
in  the  path  delay  equation.  In  this  evaluation,  the  square 
of  the  SNR  for  a  single  satellite  pass  is  used  as  the 
measurement.  The  squared  value  is  used  so  that  the 
representation  for  SNR  in  equation  (2)  can  be  used 
directly. 

A  model  for  the  antenna  gain  is  needed  in  order  to  use  the 
SNR  for  reflector  identification.  The  line  of  sight  vectors 
are  rotated  from  the  local  frame  into  the  antenna  frame 
using  an  initial  estimate  for  the  spacecraft  local-to-body 
attitude.  Given  the  line  of  sight  vectors  in  the  antenna 
frame,  the  azimuth  and  elevation  of  the  GPS  signal 
relative  to  the  antenna  is  determined,  and  thus  the  antenna 
gain  for  that  satellite  pass  can  be  found. 

The  state  vector  for  the  batch  filter  is: 

x  =  [d  az  el  A0  Am]  (4) 

The  state  is  composed  of  five  variables,  where  d  is  the 
distance  to  the  reflector,  az  and  el  are  the  azimuth  and 
elevation  of  the  plane  normal  vector,  A0  is  the  constant 
part  of  the  direct  signal  amplitude  in  amplitude  units,  and 
Am  is  the  multipath  amplitude  in  amplitude  units.  The 
azimuth  and  elevation  of  the  reflector  are  related  to  the 
plane  unit  normal  vector  by  the  following  equation: 

«  =  [sin(e/)  cos(e/)cos(oz)  cos(e/)sin(e/)],  (5) 

where  an  elevation  angle  of  90  degrees  is  defined  to  be  in 
the  positive  ^-direction  and  an  azimuth  of  zero  degrees  is 
in  the  positive  y-direction. 
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After  an  update  to  the  state  is  computed,  it  is  added  to  the 
current  estimate  of  the  state,  and  this  process  continues 
until  the  state  estimate  converges. 


Before  the  batch  estimation  algorithm  can  be  performed 
on  the  SNR  measurement  data,  the  state  variables  must  be 
initialized.  This  initialization  is  done  by  searching  over  a 
coarse  three-dimensional  grid  containing  values  of 
distance,  azimuth,  and  elevation,  representing  the  full 
range  of  possible  reflector  locations.  For  each  of  these 
assumed  locations  and  the  GPS  lines  of  sight  for  the 
satellite  pass,  a  path  delay  profile  is  computed,  using 
equation  (1).  Note  that  if  the  line  of  sight  vector  drops 
below  the  planar  reflector,  the  path  delay  is  set  to  zero  and 
no  multipath  is  considered  for  this  region.  Additionally, 
an  estimate  for  the  constant  part  of  the  direct  signal,  A0,  is 
computed  by  performing  a  least  squares  fit  to  the  SNR 
data  as  shown  in  equation  (5). 

SNR  «  A0Aa  (5) 

Once  the  path  delay  and  direct  signal  strength  are  known, 
an  estimate  of  the  multipath  amplitude,  Am,  is  computed 
by  solving  equation  (2)  using  least  squares,  where  the 
composite  amplitude  is  the  measured  SNR.  After  trying 
each  of  the  possible  reflector  locations  and  computing  the 
corresponding  direct  and  multipath  amplitudes,  the 
location  and  amplitudes  with  the  lowest  SNR  residual  is 
used  to  initialize  the  batch  algorithm. 


After  an  initial  estimate  for  the  state  is  found,  the  batch 
estimator  is  run.  The  measurement  sensitivity  is 
computed  by  taking  partials  of  the  measurement  equation, 
defined  in  equation  (2),  with  respect  to  the  state  variables. 
The  measurement  partial  matrix,  H,  is  defined  by: 


afe)  ate)  ate)  ate_)  ate)' 

dd  daz  del  dAn  dAm 


M) 


—  A„  Aa  (e  •  n)sin\|/ 


= -4  AmA0Aad(-  ev  cos  el  sin  az  +  e .  cos  el  cos  az) siny 

daz 

)=-4AmAcAad(e  x  cos  el  —  ey  sin  el  cos  az  -  e,  sin  el  sin  az)sin\|/ 


-2A0Aa +2AmAa  cosy 


2Am+2A0Aacosy 


The  measurement  residual  is  the  difference  between  the 
measured  SNR2  and  that  predicted  from  (2)  using  the 
initial  estimate  of  the  state.  The  update  to  the  state,  x ,  is 
computed  in  a  least  squares  sense  using  the  measurement 
partial  matrix  and  the  measurement  residual,  y,  as  shown 
in  equation  (7). 


HtUx  =  H  Ty 


(7) 


With  the  computed  state  vector  for  each  antenna,  carrier 
phase  corrections  are  computed  using  equation  (3)  for 
each  satellite  pass  in  which  a  reflector  is  identified.  The 
differential  phase  correction  is  computed  by  differencing 
the  corrections  for  the  individual  antennas.  The 
differential  phase  corrections  are  then  subtracted  from  the 
measured  values  and  the  updated  differential  phase 
measurements  are  filtered  using  an  attitude  estimation 
algorithm  developed  by  Ward  [Ward,  1996]. 

The  performance  of  this  algorithm  is  tested  by  comparing 
the  phase  residuals  before  and  after  the  corrections  for  an 
effective  reflector  has  been  applied.  For  flight  data,  in 
addition  to  comparing  the  residuals,  the  resulting 
improvement  in  the  attitude  solution  will  be  noted. 

RESULTS 

Simulations 

To  aid  in  the  development  of  this  approach,  multipath 
data  were  simulated  using  the  Ohio  State  Basic  Scattering 
Code.  These  simulations  were  used  to  study  multipath  in 
simple  reflecting  environments  using  ideal  antennas.  The 
data  generated  included  a  single  large  planar  reflector  and 
only  specular  reflection  was  considered.  This  was  done 
primarily  as  a  tool  to  test  the  theory  behind  the  reflector 
identification  method.  The  data  were  simulated  and  both 
SNR  and  carrier  phase  measurements  were  produced. 
The  simulated  measurements  were  shown  to  be  the  same 
as  data  calculated  using  models  presented  previously  for 
the  same  reflector  location  and  multipath  amplitude.  The 
data  were  processed  for  a  single  antenna  using  the 
reflector  identification  algorithm  and  the  reflector  and 
amplitudes  were  accurately  found. 

Ground  Tests 

In  order  to  test  the  algorithm  in  real  applications,  data 
collected  by  C.  Behre  at  the  Naval  Research  Laboratory 
test  facility  on  October  7,  1995  were  processed.  A 
Trimble  TANS  Vector  attitude  receiver  was  used  for  the 
testing  [Cohen,  1992],  [Trimble,  1993].  Three  Trimble 
patch  antennas,  a  master  and  number  1  and  3  slaves,  were 
affixed  to  a  JAWSAT  satellite  mockup  mounted  upon  a 
three-axis  rotational  spin  table.  The  antenna  height  was 
approximately  1.5  meters  above  a  concrete  surface. 
Additionally,  a  0.6  by  0.6  square  aluminum  plate 
reflector,  inclined  at  45  degrees  was  mounted  to  the 
comer  in  which  the  number  2  slave  would  normally  be 
affixed.  During  the  experiment,  90  minutes  of  static  data 
were  collected.  Upon  evaluation  the  data  were  found  to 
show  very  clear  signs  of  multipath,  both  in  the  carrier 
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phase  and  in  the  SNR.  A  photograph  of  the  experimental 
setup  is  shown  in  Figure  3. 


Figure  3:  Ground  test  experimental  setup 

A  satellite  pass  of  GPS  PRN  17,  lasting  approximately 
one  hour  was  chosen  for  evaluation.  The  reflector 
identification  algorithm  is  implemented  using  the  raw 
SNR  data  and  the  antenna  gain  data  for  the  Trimble  patch 
antenna.  Upon  implementation  of  the  batch  estimator,  it 
was  found  that  the  state  variables  corresponding  to  the 
azimuth  and  elevation  of  the  plate  are  not  observable  for  a 
single  satellite  pass  of  data.  Therefore,  the  azimuth  and 
elevation  found  during  the  initialization  processed  were 
held  fixed  during  the  batch  process,  while  the  distance  to 
the  reflector  and  the  direct  and  multipath  amplitudes  were 
estimated.  After  implementation  of  the  batch  estimator, 
the  reflector  identified  has  the  state  shown  in  Table  1. 


Table  1:  Reflector  locations  for  ground  test  data 


ant 

d  (m) 

E3H 

EHI 

A0  (AMU) 

4*  (AMU) 

m 

1.790 

WSM 

81.7 

31.059 

1.472 

i 

1.837 

9.0 

86.1 

22.964 

2.214 

3 

1.867 

92.8 

68.1 

20.927 

0.819 

Figures  4-6  show  the  measured  and  computed  SNR2  for 
the  three  antennas  for  the  states  in  Table  1.  The 
improvement  in  the  residual  phase  using  these  reflectors 
is  shown  in  Figures  7  and  8.  Upon  applying  this 
correction  to  the  phase  for  the  entire  pass,  the  RMS  of  the 
residuals  for  baseline  1  reduces  from  11.37  mm  to  8.56 
mm,  corresponding  to  a  25%  error  reduction.  The  RMS 
of  the  residuals  for  baseline  3  changes  only  slightly  from 
29.27  mm  to  29.08  mm. 


Figure  4:  Measured  (blue)  and  computed  SNR2  (green) 
for  the  master  antenna 


for  the  number  one  slave  antenna 


for  the  number  three  slave  antenna 
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phase  corrections  (green)  for  baseline  1 


phase  corrections  (green)  for  baseline  3 

For  these  ground  test  data,  the  method  works  very  well  in 
identifying  an  effective  reflector  by  which  to  model 
multipath.  The  ground  was  identified  as  the  primary 
source  of  multipath.  The  sources  of  error  for  the  results 
shown  in  Figures  7  and  8  are  believed  to  be  caused  by  two 
things.  First,  the  multipath  amplitude  is  being  modeled  as 
a  constant  value.  An  adaptive  means  by  which  to 
estimate  the  multipath  amplitude  would  improve  these 
results,  since  it  is  clear  in  both  the  SNR  and  phase  plots 
that  the  multipath  amplitude  does  vary.  Additionally, 
there  is  a  lower  frequency  oscillation  in  the  master  and  the 
number  3  slave  antennas  that  is  not  being  accurately 
modeled  by  the  single  large  reflector  approximately  1.75 
meters  away.  Estimation  of  a  second  effective  reflector, 
which  is  possibly  the  square  aluminum  plate,  would  likely 
also  decrease  the  differential  phase  residuals. 

CRISTA-SPAS 

The  Space  Shuttle  Atlantis  deployed  the  CRISTA-SPAS 
satellite  on  November  3,  1994.  A  photograph  of  the 
satellite  being  released  is  shown  in  Figure  9.  During  its 
11 -day  mission,  it  remained  within  a  few  kilometers  of 
the  Shuttle,  with  the  primary  purpose  of  the  mission  being 


observation  of  the  Earth’s  middle  atmosphere.  Three-axis 
control  of  the  spacecraft  was  required  in  order  to  point  a 
telescope  at  an  area  62.9  kilometers  above  the  WGS-84 
ellipsoid.  A  star  tracker-gyro  inertial  reference  unit  (IRU) 
provided  attitude  information  for  closed  loop  control  and 
an  Alcatel  GPS  receiver  provided  position  information 
[Brock,  et  al,  1995]. 


Figure  9:  Photograph  of  CRISTA-SPAS  (Courtesy  J. 
Rodden,  Loral  Space  Systems) 


Before  the  CRISTA-SPAS  data  are  processed  using  the 
reflector  identification  method,  the  spatial  repeatability  of 
the  data  is  tested.  This  is  done  by  identifying  four  GPS 
tracks  that  follow  the  same  path.  The  azimuth  and 
elevation  of  the  four  tracks  are  shown  in  Figure  10.  The 
SNR  and  residual  differential  phase  for  these  satellite 
passes  are  shown  in  Figures  1 1  and  12,  respectively. 
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Figure  10:  4  satellite  passes  following  similar  tracks 


Figure  1 1 :  SNR  for  the  Slave  3  antenna  for  the  4  passes 


Figure  12:  Residual  differential  phase  for  baseline  1  for 
the  4  passes 

Since  both  the  SNR  and  residual  phase  follow  the  same 
trends  for  the  similar  satellite  passes,  it  can  be  inferred 
that  the  data  are  highly  spatially  correlated.  It  is 
reasonable,  therefore,  to  assume  that  the  reflector 
identification  method  will  be  effective  in  isolating  an 
effective  reflector  location  that  will  reduce  the  phase 
residuals  for  not  only  a  single  pass,  but  other  passes  that 
follow  the  same  path  in  the  sky.  Also,  once  an  effective 
reflector  has  been  found  for  that  portion  of  the  sky,  there 
is  no  need  to  re-solve  for  a  new  location.  The  previous 
computed  residuals  may  simply  be  applied  or  the  reflector 


location  may  be  iterated  upon  by  using  the  previous 
reflector  location  as  the  a  priori  state. 

Another  test  of  the  repeatability  of  phase  multipath  was 
conducted  by  implementing  a  phase  map  correction.  A 
detailed  explanation  of  the  process  and  results  are  shown 
in  [Axelrad  and  Reichert,  1997].  This  involves 
computing  phase  residuals  for  all  32  hours  of  satellite  data 
and  generating  a  phase  correction  map  as  a  function  of  the 
line  of  sight  vector.  This  is  based  on  the  method 
presented  in  [Cohen  and  Parkinson,  1991].  In  this 
evaluation,  a  number  of  methods  were  used  to  implement 
a  fit  to  the  data,  including  a  spherical  harmonic  fit,  a 
polynomial  fit,  and  a  correction  grid  method.  The 
correction  grid  fit  was  shown  to  provide  the  best  fit  to  the 
phase  residuals. 

The  correction  grid  is  generated  by  sorting  the  data  into 
bins  based  on  azimuth  and  elevation  values  of  the 
incoming  signal.  The  correction  value  for  each  grid  point 
is  set  to  the  mean  value  of  all  the  residuals  within  each 
bin.  The  grid  computed  has  a  rectangular  resolution  of 
approximately  1  deg  x  1  deg  and  is  shown  for  baseline 
one  in  Figure  13.  The  correction  grid  may  either  be 
applied  as  either  a  table  look-up  or  it  may  be  inteipolated. 
Upon  implementation  of  the  correction  grid  to  the  phase 
data  using  interpolation,  the  RMS  of  the  residuals  were 
reduced  by  51%,  52%,  and  41%  for  baselines  1,  2,  and  3, 
respectively.  A  linearized  Kalman  filter  is  implemented 
using  the  improved  phase  measurements.  The  resulting 
attitude  error  improves  from  0.12  to  0.05  degrees  in  yaw, 
0.10  to  0.04  degrees  in  roll  and  0.20  to  0.14  degrees  in 
pitch. 


1  Degree  G ridded  Residual  Phase  (mm)  for  Baseline  1 


Figure  13:  Phase  correction  grid  map  for  Baseline  1 


This  method  of  computing  corrections  to  the  phase,  using 
the  attitude  solution  is  somewhat  flawed  because  it 
requires  ideal  attitude  and  baseline  solutions  in  order  to 
accurately  determine  the  phase  residuals.  Any  error  in  the 
attitude  or  baselines  directly  maps  into  the  computation 
for  the  residual  phase.  Therefore  it  is  difficult  to  separate 
errors  in  attitude  from  multipath.  In  order  to  avoid  the 
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problem  a  dependence  on  SNR,  rather  than  phase,  is 
utilized  because  the  SNR  is  much  less  sensitive  to  attitude 
errors. 

Utilizing  the  SNR,  the  reflector  identification  algorithm  is 
implemented  for  the  CRISTA-SPAS  data.  The  gain 
pattern  used  for  the  ground  tests  for  the  Trimble  patch 
antennas  did  not  match  well  with  the  received  SNR 
values.  Therefore  a  simple  two-dimensional  polynomial 
fit  to  the  SNR  data  was  used  for  the  antenna  gain  value, 

A* 

The  reflector  identification  method  was  implemented  for 
a  number  of  satellite  passes.  The  SNR  fits  for  the  four 
antennas  for  a  single  pass  of  satellite  PRN  1  are  shown  in 
Figures  14  -  17.  One  successful  pass  is  for  the  off- 
baseline  formed  by  the  number  2  and  3  slave  antennas. 
The  phase  correction  for  this  off-baseline  is  shown  in 
Figure  18.  Unfortunately,  an  appropriate  effective 
reflector  was  not  located  for  the  master  antenna,  as  is 
represented  by  the  phase  residual  for  baseline  1  shown  in 
Figure  19,  so  the  phase  measurements  could  not  be 
adequately  updated  for  this  satellite  pass. 


fit  for  the  Master  antenna 


fit  for  the  Slave  1  antenna 


fit  for  the  Slave  2  antenna 


fit  for  the  Slave  3  antenna 


baseline 
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Tables  2  and  3  list  results  for  three  baselines  for  PRN  1, 
baseline  2  and  PRN  25,  baselines  1  and  2,  1  in  which  the 
residuals  are  improved.  Table  2  includes  the  reflector 
location  and  the  multipath  amplitude  and  Table  3  shows 
the  residual  differential  phase  before  and  after  the  phase 
corrections  are  applied.  For  these  satellite^aseline 
combinations  in  which  phase  residuals  were  improved, 
the  phases  measurements  were  corrected  and  an  attitude 
point  solution  was  recomputed  using  the  algorithm 
developed  by  Ward  [Ward,  1996].  Because  only  a  total  of 
three  baselines  were  corrected  with  minimal  phase 
improvements,  no  appreciable  attitude  improvement  was 
found. 


Table  2:  Reflector  location  and  amplitude  for 
CRISTA-SPAS  data _ 


PRN 

antenna 

d 

(m)  ... 

az 

(deg) 

el 

(deg) 

d;n 

(AMU) 

1 

master 

.078 

285 

75 

2.45 

1 

slave  2 

.079 

240 

72 

3.84 

25 

master 

.428 

220 

80 

0.69 

25 

slave  1 

.385 

300 

70 

1.87 

25 

slave  2 

.141 

80 

70 

1.31 

Table  3:  Phase  improvements  for  CRISTA-SPAS  data 


PRN 

Baseline 

residual  before 

residual  after 

1 

m-2 

7.5095 

7.3611 

25 

m-1 

6.4085 

5.5335 

25 

m-2 

6.7899 

6.4252 

CONCLUSIONS  AND  FUTURE  WORK 


In  this  paper  a  method  to  reduce  carrier  phase  multipath 
by  identifying  an  effective  reflector  using  the  GPS  SNR 
measurement  was  presented.  The  method  utilizes  the 
known  spatial  correlation  of  multipath  to  identify  a 
reflector  location  and  direct  signal  and  multipath 
amplitudes. 

The  results  shown  for  the  ground  data  indicate  significant 
improvement  to  the  residual  differential  phase  using  only 
a  single  reflector  and  multipath  amplitude  identified  from 


a  batch  estimation  scheme.  It  is  likely  that  more 
improvement  would  result  if  the  multipath  amplitude 
were  to  be  adaptively  estimated  over  time  to  better 
account  for  the  changing  reflectivity  of  the  reflecting 
surface.  The  results  found  by  implementing  the  method 
for  the  CRISTA-SPAS  flight  data  indicate  some  potential 
to  identify  a  reflector  capable  of  reducing  carrier  phase 
multipath.  However,  these  data  demonstrated  that  the 
ability  to  identify  a  reflector  from  SNR  data  is  limited  by 
the  poor  knowledge  of  the  antenna  gain  pattern. 

In  future  work  the  batch  filter  will  be  replaced  by  an 
adaptive  filter  which  will  allow  the  reflector  location  and 
multipath  amplitude  to  slowly  change  as  the  satellite 
passes  overhead.  Additionally,  a  second  reflector  will  be 
identified  in  some  situations  where  it  is  clear  that  both 
high  and  low  frequency  multipath  components  exist, 
being  caused  by  reflectors  both  far  from  and  near  to  the 
antenna. 
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Chapter  1 


Introduction 


Accurate  relative  navigation  of  spacecraft  permits  capabilities  such  as  rendezvous  and  formation 
flying — key  aspects  of  many  current  and  future  space  missions.  In  low  earth  orbit  (LEO),  such 
missions  can  use  the  Global  Positioning  System  (GPS)  signal  for  relative  positioning  and  data  time¬ 
tagging.  However,  little  has  been  published  on  the  natural  extension  to  GPS  relative  navigation: 
precise  time  and  frequency  transfer  on  orbit.  In  terrestrial  time  keeping,  GPS  is  widely  used 
for  nanosecond-  and  even  picosecond-level  time  transfer,  given  long  averaging  times  and  accurate 
knowledge  of  positions  and  GPS  signal  error  sources.  My  proposed  research  will  investigate  methods 
for  extending  GPS  time  transfer  to  satellite-to-satellite  applications  at  the  five  nanosecond  level 
or  better  in  near  real  time.  The  primary  topic  is  determination  and  demonstration  of  algorithms 
for  precise  relative  position,  velocity,  and  GPS  receiver  clock  bias  and  frequency  through  use  of 
both  GPS  code  and  carrier  phase  measurements.  The  algorithms  will  be  tested  using  flight  data 
and  hardware-in-the-loop  simulation  data.  This  report  discusses  the  progress  to  date  and  the  steps 
proposed  to  reach  the  stated  goals. 
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CHAPTER  1.  INTRODUCTION 


1.1  Previous  Research 

1.1.1  GPS  Relative  Navigation 

Recent  years  have  seen  an  abundance  of  GPS  relative  navigation  research.  An  extensive  evalu¬ 
ation  of  satellite-to-satellite  GPS  relative  navigation  algorithms  is  presented  by  Binning  [1].  He 
generates  data  for  two  co-orbiting  satellites  in  LEO  using  both  software  and  hardware-in-the-loop 
simulations  and  presents  the  results  for  four  separate  algorithms.  Two  of  his  methods  compute  the 
relative  state  by  subtracting  absolute  orbit  estimates  determined  using  GPS  P-code  pseudorange 
observations,  while  the  two  other  methods  explicitly  estimate  the  relative  state.  The  first  method 
generates  independent  orbit  estimates  of  each  vehicle  using  all  the  GPS  satellites  in  view  to  the 
respective  receivers,  achieving  a  relative  position  accuracy  of  about  1  m.  A  more  complex  method 
interpolates  the  pseudorange  observations  from  the  two  receivers  to  a  common  time  and  solves  for 
the  absolute  orbits  using  only  common-view  satellites.  An  extension  to  this  method  uses  single 
difference  pseudoranges  to  estimate  the  relative  state  directly,  while  also  estimating  the  absolute 
state  of  one  of  the  vehicles  as  before.  These  two  common-view  techniques  achieve  relative  position 
accuracies  of  approximately  15-25  cm.  Finally,  Binning  presents  an  algorithm  which  uses  both 
GPS  transmission  frequencies  to  create  a  widelane  carrier  phase  observation.  Using  single  differ¬ 
ences  with  this  type  of  observation  increases  accuracy  to  approximately  3  cm  in  relative  position. 
It  should  be  noted  that  Binning  did  not  use  a  more  accurate  external  reference  clock  during  the 
hardware  simulations,  so  the  accuracy  of  his  clock  bias  and  frequency  estimates  are  not  known. 

Additional  simulation  work  has  been  performed  by  Upadhyay  et  al.  [2]  in  support  of  automated 
rendezvous  and  docking  studies  for  NASA.  Their  two-stage  GPS  relative  navigation  scheme  con¬ 
sists  of  a  meter-level  filter  and  a  centimeter-level  filter.  The  meter-level  filter  generates  smoothed 
absolute  and  relative  orbit  estimates  for  two  vehicles  using  GPS  pseudorange  and  either  Doppler 
or  range-rate.  The  centimeter-level  filter  refines  the  relative  state  estimates  by  solving  for  the 
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double-differenced  carrier  phase  ambiguities.  Hardware-in-the-loop  tests  for  the  meter-level  filter 
produces  RSS  relative  position  and  velocity  errors  of  0.887  m  and  0.8  cm/s,  respectively.  Results 
from  software  simulated  data  on  the  centimeter- level  filter  show  relative  RSS  accuracies  of  0.11  m 
in  position  and  1  cm/s  in  velocity. 

Cox  and  Brading  [3]  also  developed  a  GPS  relative  navigation  filter  that  resolves  the  carrier 
phase  ambiguities.  In  their  paper,  they  discuss  implementation  of  the  LAMBDA  method  for  ambi¬ 
guity  resolution  in  an  optimal  Kalman  filter  and  evaluate  its  performance  on  simulated  data.  The 
measurement  types  used  are  single-frequency,  uncorrupted  P-code  and  double-differenced  accumu¬ 
lated  phase.  For  the  10-km  baseline  case,  the  relative  radial  position  errors  are  generally  less  than 
2  cm.  Increasing  the  baseline  to  100  km  results  in  radial  position  errors  of  approximately  3-10  cm. 

Two  notable  flight  experiments  utilizing  GPS  for  relative  navigation  occurred  with  the  NASA 
STS-80  Shuttle  mission  and  the  NASDA  ETS-VII  flight  experiment.  Schiesser  et  al.  [4]  present 
results  from  the  Shuttle/ORFEUS-SPAS  free-flyer  rendezvous  during  the  STS-80  mission  in  De¬ 
cember  1996.  Their  filter  uses  GPS  coarse/acquisition  (C/A)  code  pseudoranges  to  solve  for  the 
absolute  orbits  of  the  two  vehicles,  then  differences  them  to  form  the  relative  orbit.  The  filter  is  de¬ 
signed  to  processes  asynchronous  observations  from  the  two  receivers  in  real  time,  accepting  input 
from  the  Shuttle  IMU  for  maneuver  compensation.  A  unique  element  of  this  filter  is  that  it  uses 
the  same  filter-estimated  range  bias  state  and  covariance  matrix  elements  for  both  vehicles  when 
processing  observations  from  common-view  GPS  satellites.  Given  the  highly  dynamic  environment 
during  rendezvous,  the  filter  obtained  RSS  relative  position  and  velocity  errors  of  9  m  and  12.6 
cm/s,  respectively. 

Most  recently,  D ’Souza  et  al.  [5]  from  Draper  Lab  developed  a  filter  for  validation  of  algorithms 
for  use  with  two  NASDA  projects,  Engineering  Test  Satellite  7  (ETS-VII)  and  the  H-II  Transfer 
Vehicle  (HTV).  The  filter  estimates  explicitly  the  relative  position,  velocity,  and  GPS  receiver  clock 
bias  and  frequency  using  single  difference  C/A  code  pseudoranges.  The  asynchronous  pseudoranges 
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from  one  receiver  are  linearly  interpolated  to  the  times  of  the  other  to  allow  measurement  differ¬ 
encing.  Excellent  simulation  results  were  obtained.  Presuming  NASDA  implemented  this  filter  for 
ETS-VII,  Mokuno  et  al.  [6]  indicate  that  on-orbit  results  were  also  good,  obtaining  a  GPS-derived 
relative  range  within  1  m  of  the  rendezvous  radar  measurement  and  an  in-plane  velocity  error 
within  2  cm/s. 

As  for  GPS  relative  navigation  combined  with  time  transfer  for  orbiting  platforms,  very  few 
references  appear  in  the  literature.  Cox  and  Brading  have  published  two  such  papers  [7,  8].  In  the 
most  recent  of  the  two  (1997),  they  use  software  simulated  orbits  with  baselines  of  100  km  and  10 
km  separating  the  two  vehicles.  (The  same  simulation  is  used  for  their  subsequent  paper  using  the 
LAMBDA  method  [3],  but  there  is  no  mention  of  time  transfer  in  that  paper.)  Dual-frequency, 
uncorrupted  P(Y)  code  and  carrier  measurements  are  provided  by  the  simulation.  The  filter  creates 
differences  and  sums  of  pseudoranges  and  differences  of  accumulated  phase  measurements  to  use 
as  filter  observations.  The  filter  estimates  absolute  position,  velocity,  clock  time  error,  and  clock 
frequency  error  for  each  vehicle  (16  states),  as  well  as  floating-point  single-difference  carrier  phase 
ambiguities  for  each  common-view  GPS  satellite.  At  100  km,  the  relative  RMS  accuracies  are 
given  as  7  cm  in  position,  0.01  cm/s  in  velocity,  0.4  ns  in  time  (ignoring  common-mode  errors), 
and  7  x  10-12  in  fractional  frequency.  Note,  however,  that  the  simulation  employs  a  second-order 
stochastic  model  to  generate  the  true  receiver  clock  errors,  and  that  the  second-order  clock  model  in 
the  filter  uses  the  same  parameters  as  the  truth  for  propagating  clock  errors.  Hardware-in-the-loop 
simulations  do  not  have  the  benefit  of  setting  or  knowing  precisely  the  clock  behavior. 


1.1.2  GPS  Time  Transfer 


Although  not  much  has  been  published  on  GPS  time  transfer  in  orbit,  extensive  research  has  been 
performed  on  terrestrial  applications  of  GPS  time  transfer.  Indeed,  use  of  GPS  is  currently  the 
foremost  method  for  precise  time  and  time  interval  distribution.  Detailed  summaries  of  precise 
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timing  using  GPS  can  be  found  in  [9,  10,  11].  Allan  and  Weiss  first  proposed  the  now  standard 
GPS  common- view  (CV)  technique,  which  involves  differencing  C/A  code  pseudorange  observations 
from  different  timing  sites  with  known  locations  [12].  Current  multi-channel  CV  results  by  Levine 
indicate  a  3-5  ns  RMS  accuracy  [13].  Larson  and  Levine  further  refine  time  transfer  accuracies  by 
employing  the  GPS  carrier  phase,  achieving  a  stability  of  100  ps  and  frequency  uncertainty  of  2 
parts  in  1015  for  1  day  averages  [14]. 

On-orbit  time  transfer  will  not  have  the  benefit  of  extremely  precise  knowledge  of  the  antenna 
locations,  as  in  most  terrestrial  time  transfer  applications.  Error  sources  must  be  treated  differently 
in  orbit,  and  the  averaging  times  must  be  shorter  for  meaningful  application  in  real  time.  However, 
picosecond  accuracies  are  not  required  for  the  on-orbit  application.  Thus,  a  combination  of  precise 
relative  navigation  and  time  transfer  techniques  is  sought. 


1.2  Contributions 

The  primary  contribution  of  this  research  will  be  the  ability  to  synchronize  the  clocks  of  a  pair  of 
orbiting  satellites  in  real  time  to  five  nanoseconds  with  respect  to  each  other.  Previous  research  has 
shown  precise  relative  navigation  of  satellites  and,  separately,  precise  time  and  frequency  transfer 
using  GPS.  However,  none  has  shown  precise  relative  navigation  and  timing  of  satellites  in  orbit  in 
real  time  using  actual  hardware  data. 

Another  contribution  of  this  research  is  the  development  of  a  time  transfer  test  bed  for  per¬ 
forming  hardware  simulations  of  multiple  satellites.  This  research  is  being  performed  at  the  Naval 
Research  Laboratory,  which  is  uniquely  equipped  for  such  a  test  bed.  Details  given  in  a  later  sec¬ 
tion  show  that  the  test  bed  will  measure  the  absolute  and  relative  clock  behavior  of  the  hardware, 
providing  a  truth  for  comparison.  Thus,  the  relative  timing  algorithms  can  be  designed  to  reach 
an  estimation  accuracy  equivalent  to  the  measured  precision  of  the  hardware. 
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1.3  Overview  of  Report 

Chapter  2  provides  the  details  of  the  current  version  of  the  extended  Kalman  filter  for  relative 
navigation.  Chapter  3  describes  the  flight  experiment  performed  during  a  Shuttle  mission  and  the 
results  obtained  with  the  current  filter.  Chapter  4  shows  the  performance  of  the  filter  on  a  data 
set  which  more  closely  resembles  expected  spacecraft  dynamics  and  data  quality.  These  data  were 
obtained  from  a  previous  hardware  simulation  run  by  Binning  in  [1].  Finally,  Chapter  5  details  the 
proposed  time  transfer  test  bed  and  testing  procedure. 


Chapter  2 


Relative  State  Kalman  Filter 


2.1  Overview 

An  extended  Kalman  filter  that  directly  estimates  the  relative  state  between  two  co-orbiting  ve¬ 
hicles  using  single  difference  pseudorange  observations  has  been  developed.  This  filter  is  used  in 
conjunction  with  an  existing  GPS  absolute  orbit  determination  utility,  which  provides  an  estimate 
of  the  absolute  orbit  of  the  passive  vehicle  to  the  relative  filter.  Use  of  single  difference  pseudoranges 
reduces  errors  inherent  in  using  GPS  by  removing  common  mode  errors,  such  as  the  satellite  clock 
error.  Furthermore,  directly  estimating  a  relative  state  allows  the  use  of  the  Clohessy- Wiltshire 
(CW)  equations  for  state  and  error  propagation.  In  fact,  in  its  current  form  the  relative  filter 
requires  no  integration. 

Specifically,  this  filter  estimates  a  relative  position  and  velocity  and  relative  receiver  clock  bias 
and  frequency.  It  uses  the  CW  equations  with  velocity  aiding  to  propagate  the  position  and  velocity, 
incorporating  the  absolute  orbit  estimate  of  the  passive  vehicle  to  compute  the  mean  motion.  The 
CW  equations  are  also  used  to  propagate  the  covariance.  The  measurement  update  computes  the 
observed  and  predicted  single  difference  pseudoranges,  where  the  predicted  quantity  accounts  for 
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the  body-fixed  GPS  antenna  offset  from  the  center  of  gravity  by  using  a  reported  attitude  estimate. 
The  relative  state  sr  is  expressed  by 


Sr  — 


br 


(2.1) 


where  rr  is  position,  vr  is  velocity,  br  is  clock  bias,  and  fr  is  frequency.  The  filter  position  and 
velocity  coordinate  frame  is  centered  on  one  of  the  vehicles  and  expressed  in  radial,  in-track,  and 
cross-track  (RIC)  directions.  The  relative  state  is  referenced  to  the  passive  vehicle  in  the  case  of  a 
rendezvous.  In  the  case  of  two  non-thrusting  vehicles,  one  is  arbitrarily  chosen  to  define  the  frame 
of  reference.  The  covariance  and  process  noise  values  are  also  specified  in  the  RIC  directions.  This 
formulation  permits  more  direct  use  of  the  CW  equations  for  state  and  covariance  propagation. 

The  Relative  State  Kalman  Filter  (RSKF)  is  designed  to  handle  asynchronous  GPS  code  pseu¬ 
doranges  from  the  two  receivers,  absolute  orbit  estimates  for  the  passive  vehicle,  and  GPS  satellite 
position  vectors.  Since  the  RSKF  is  designed  to  work  in  tandem  with  an  external  GPS  absolute 
orbit  determination  (AOD)  utility,  output  from  it  is  used  as  input  to  the  RSKF.  The  external  AOD 
program  generates  the  absolute  orbit  estimates  of  the  passive  vehicle  using  GPS  code  pseudoranges 
and  broadcast  GPS  satellite  orbits.  The  orbit  estimates,  along  with  the  GPS  observations  and 
satellite  positions  used  to  compute  them,  are  output  from  the  AOD  program  for  use  as  input  to 
the  RSKF.  Within  the  RSKF  the  active  vehicle  pseudoranges  are  interpolated  to  the  times  of  the 
passive  vehicle  pseudoranges  and  positions. 

Inertial-to-body  attitude  quaternions  and  delta-V  vectors  can  also  be  input  to  the  RSKF.  The 
attitude  quaternions  are  needed  to  account  for  the  GPS  antenna  CG  offsets,  which  are  specified  in 
a  body-fixed  frame.  The  inertial-to-RIC  rotation  is  computed  using  the  absolute  position  estimate 
of  the  passive  vehicle.  The  delta-V  vectors  are  used  in  the  time  update  propagation  to  account  for 
maneuvers  above  a  specified  threshold. 
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2.2  Time  Update 


As  stated  previously,  the  CW  equations  are  used  in  this  filter  for  relative  state  and  covariance 
propagation.  From  Prussing  and  Conway  [15],  the  CW  equations  in  state  transition  matrix  form 
are 

4-3  c  00  f  0 

6(s-nAt)  1  0  4s=2nAt  0 

0  0  c  0  Of 

3ns  0  0  c  2s  0 

— 6n(l  -  c)  0  0  -2s  4c  -  3  0 

0  0  —ns  0  0  c 

where  n  signifies  the  mean  motion  of  the  passive  vehicle,  At  is  the  time  interval  for  propagation, 
and  s  =  sin(nAi)  and  c  =  cos(nAt).  In  the  filter,  the  mean  motion  in  the  equation  above  is  actually 
computed  as  the  angular  rotation  rate  (0)  at  the  current  time.  In  other  words,  the  mean  motion 
is  chosen  at  each  epoch  to  be  the  current  along-track  velocity  of  the  passive  vehicle  divided  by  its 
radius.  This  provides  a  more  accurate  propagation  for  elliptical  orbits  as  long  as  the  measurement 
interval  is  not  large. 


*P,v(Ai)  = 


The  relative  clock  frequency  is  assumed  constant  over  the  time  interval,  so  the  clock  terms  are 


propagated  with  [16] 

1  At 

*6,/(A<)  = 

0  1 

Thus,  the  full  relative  state  transition  matrix  has  the  form 


(2.3) 


$  (At)  = 


$r,v  (At)  06x2 

02x6  $(>,/  (At) 


(2.4) 


Incorporated  into  the  state  propagation  is  velocity  aiding,  accomplished  by  using  the  Av  vectors 
of  the  active  vehicle.  The  velocity  vectors  are  added  directly  into  the  linear  state  propagation  at  the 
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given  time  as  impulsive  maneuvers,  with  potentially  more  than  one  maneuver  between  measurement 
epochs.  The  following  example  illustrates  the  equations  for  the  occurrence  of  one  Av  during  the 
interval.  As  such,  the  state  must  be  propagated  from  4  to  tk+i  with  a  maneuver  at  tAv  such  that 
tk  <  tAv  <  tk+ 1.  In  this  case,  the  state  vector  is  propagated  as  follows: 


sAv  =  $(4^,4)  s£ 

(2.5) 

SA«  =  sAv  +  [°3xl  Av-  OLlf 

(2.6) 

Sfc+1  =  ^  (4+l>  ^Au)  sAu 

(2.7) 

where  the  ~  and  +  indicate  state  vector  values  before  or  after  the  measurement  update  or  a  ma¬ 
neuver.  The  distinction  is  clear  from  the  context. 

Not  portrayed  in  the  above  equations  is  the  fact  that  the  CW  state  propagation  is  performed 
in  cylindrical  coordinates,  with  the  along-track  direction  being  along  the  arc  of  the  orbit.  Since 
this  filter  operates  in  the  rectangular  RIC  coordinates,  a  rotation  is  performed  to  cylindrical  RIC 
coordinates  immediately  before  each  state  propagation  step  (e.g.,  Eq.  2.5),  then  back  to  rectangular 
immediately  after.  The  equations  for  such  a  conversion  are  given  in  [17].  However,  for  this  filter, 
those  equations  are  re-derived  for  the  case  that  r  is  not  zero,  providing  a  slight  increase  in  accuracy. 

Ignoring  any  potential  maneuvers,  the  covariance  is  propagated  over  the  entire  measurement 
interval  using 

P*+1  =  *  (4+i ,4)  Pfc  *  (4+i ,  tkf  +  Qk  (2-8) 

The  process  noise  matrix  Q  is  derived  from  Gelb,  equation  3.6-13  [18], 

tfc+i 

TfcQfclf  =  /  (4+i ) T)  G  (r)  Q  (t)  Gt  (t )  3>t  (4+1  > T)  dr  (2.9) 

tk 

where  Q (t)  is  defined  to  be  the  spectral  density  matrix,  and  the  corresponding  discrete  state 
equation  is 


Sfc+i  =  3»fcSk  +Tkwk. 


(2.10) 
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However,  since  the  CW  equations  are  used  for  state  propagation,  w*  is  assumed  to  be  zero  mean, 
and  rfc  is  assumed  to  be  identity.  If  equal  weighting  is  assumed,  G (t)  becomes  an  identity  matrix, 
and  Eq.  2.9  reduces  to 

tk+l 

Q  k=  J  $(tk+l,T)Q(T)<f>T  (tk+UT)dT.  (2.11) 

h 

Eq.  2.11  is  solved  symbolically  using  Eq.  2.4  and  a  diagonal  representation  of  the  spectral  density 
matrix  Q  (t). 


2.3  Measurement  Update 

The  RSKF  performs  a  measurement  update  for  each  GPS  observation  separately,  so  the  equations  in 
this  section  refer  to  quantities  for  one  common  view  GPS  satellite.  Given  pseudorange  observations 
to  the  same  GPS  satellite  from  the  active  receiver  pA  and  passive  receiver  pp,  the  single  difference 
observation  is  simply 

Psd  =  Pa~  Pp-  (2.12) 

The  predicted  single  difference  pseudorange  psd  is  computed  from  the  geometric  relative  range, 
relative  clock  bias  estimate,  and  antenna  CG  offsets.  Because  the  pseudoranges  in  the  above 
equation  are  measured  from  the  GPS  antennas,  psd  is  not  necessarily  referenced  to  the  centers 
of  gravity  of  the  two  vehicles.  Therefore,  the  prediction  of  the  single  difference  pseudorange  must 
account  for  the  antenna  CG  offset  vector  rANT.  In  the  equations  below,  rANT  can  be  viewed  as 
the  sum  of  the  antenna  offset  vectors  from  the  two  vehicles.  This  assumes  similar  attitudes  for  the 
vehicles  to  promote  common  GPS  constellation  visibility. 

The  geometric  relative  range  is  found  by  differencing  the  ranges  to  the  active  and  passive  vehicles 
using  the  GPS  satellite  position  Rgps ,  passive  vehicle  absolute  position  estimate  r p,  and  predicted 
relative  position  estimate  trel  from  the  time  update.  The  computed  ranges  to  the  active  and 
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passive  vehicle  receivers  are 


fP  =  |Rgps  -  rp| 
ta  =  |Rgps  -  (rP  +  *rel  +  i/uvr)!  • 

The  predicted  single  difference  is  then  computed  from 

PSD  =  rA~  rp  +  bREL- 

The  single  difference  pseudorange  residual  is  formed  using 

Z  =  PSD  ~  PSD- 


(2.13) 

(2.14) 


(2.15) 


(2.16) 


The  partial  differential  of  Eq.  2.15  with  respect  to  the  relative  state  (Eq.  2.1)  produces  the 
measurement  gradient  matrix  H  with  the  first  three  elements  being 


Hn  = 
#12  = 
#13  = 


Xgps  —  (xp  +  xrel  +  xant) 


r  A 


YgPS  -  ( yp  +  VREL  +  VANT ) 


rA 


ZgPS  -  {zp  +  Zrel  +  Zant) 


TA 


(2.17) 

(2.18) 
(2.19) 


so  that  the  full  matrix  is 


H  = 


#11  #12  #13  0  0  0  1  0 


The  Kalman  gain  is  computed  from  the  equation 

K  =  Pl+1Hr  (HP»+lHT  +  R )"' 


(2.20) 


(2.21) 


and  the  covariance  is  updated  with 

P++1  =  (I  -  KH)  P^+1  (I  -  KHf  +  KRKt.  (2.22) 


Since  the  measurement  update  is  performed  for  each  GPS  satellite  separately,  the  measurement 
noise  matrix  ft,  is  1  x  1,  so  the  inversion  in  Eq.  2.21  is  a  simple  scalar  inversion. 
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Finally,  the  relative  state  is  updated  by 

4+i  =  sk+i+K  z-  (2-23) 

The  measurement  update  is  repeated  for  each  common  view  GPS  satellite.  After  exhausting  all  the 
observations  for  this  epoch,  the  next  measurements  are  chosen  and  the  entire  process  is  repeated. 
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Chapter  3 


Flight  Test  Results 


3.1  STS-80  Flight  Experiment 


The  STS-80  mission  launched  November  19,  1996,  and  returned  to  Earth  December  7,  1996.  During 
the  mission,  the  Space  Shuttle  Columbia  released  the  free-flying  Orbiting  and  Retrievable  Far  and 
Extreme  Ultraviolet  Spectrograph  Shuttle  Pallet  Satellite  (ORFEUS  SPAS).  Both  the  SPAS  and 
the  Shuttle  had  GPS  receivers  on  board,  which  were  used  to  record  concurrent  observations  for 
demonstration  of  relative  navigation  using  GPS.  The  SPAS  had  a  nine-channel  Laben  Tensor  GPS 
receiver  with  four  antennas,  while  the  Shuttle  used  a  six-channel  Trimble  Advanced  Navigation 
System  Quadrex  GPS  receiver  with  two  antennas  mounted  across  the  payload  bay  [4]. 

The  data  set  used  to  test  the  RSKF  comprises  the  final  73  minutes  of  rendezvous,  up  to  the 
capture  (grapple)  of  the  SPAS  with  the  Shuttle  robotic  arm.  Several  data  sets  were  taken  through¬ 
out  the  mission,  but  only  the  rendezvous  portion  had  a  valid  relative  orbit  truth  reference.  This 
segment  of  the  experiment  employed  a  laser  on  the  Shuttle  for  determining  range  and  range-rate  to 
the  SPAS.  These  data  were  converted  into  a  best  estimated  relative  trajectory  (BERT)  in  radial, 
in-track,  and  cross-track  (RIC)  coordinates. 
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Because  of  the  frequent  thrusting  of  the  Shuttle  during  rendezvous,  the  RSKF  delta-V  compen¬ 
sation  is  employed.  The  delta-V  measurements  were  produced  by  the  Shuttle  IMU  accelerometer, 
which  reported  accumulated  sensed  velocity  [4].  Given  Shuttle  delta-V  vectors  expressed  in  RIC 
coordinates,  implementation  of  velocity  aiding  in  the  filter  is  straightforward. 

The  rendezvous  phase  of  the  flight  was  performed  in  such  as  way  as  to  maximize  the  potential 
for  common  view  satellites.  The  Shuttle  payload  bay  remained  pointed  at  the  SPAS  so  the  crew 
could  maintain  visual  contact;  consequently,  the  SPAS  was  commanded  to  pitch  at  a  constant  rate 
so  that  its  GPS  antennas  pointed  in  approximately  the  same  direction  as  those  on  the  Shuttle 
throughout  the  rendezvous  [4]. 


3.2  Measurements 

Both  GPS  receivers  reported  C/A  code  phase,  LI  carrier  phase,  and  LI  Doppler.  The  RSKF  uses 
the  C/A  code  pseudoranges  from  each  common- view  GPS  satellite  to  create  the  single  difference 
observation.  The  receivers  were  not  synchronized,  so  the  measurement  time  tags  are  different. 
Therefore,  the  pseudoranges  from  the  active  vehicle  (Shuttle)  receiver  are  interpolated  to  the  times 
of  the  passive  vehicle  (SPAS)  pseudoranges  using  polint.c,  a  polynomial  interpolation  routine  from 
Numerical  Recipes  [19].  Polint.c  uses  Neville’s  algorithm  to  construct  the  interpolating  polynomial. 
As  implemented  in  the  filter,  the  desired  point  is  computed  between  the  second  and  third  values  of 
a  four-point  interpolation. 

It  should  be  noted  that  both  receivers  maintain  a  clock  bias  between  -0.5  and  +0.5  ms,  so  the 
relative  clock  behavior  for  these  receivers  has  steps  of  approximately  1  ms.  The  RSKF  searches  for 
these  jumps  by  comparing  the  current  single  difference  pseudorange  to  the  previous  one,  which  axe 
approximately  300  km  different  if  a  receiver  clock  correction  occurred.  Also,  neither  the  receivers 
nor  the  RSKF  corrects  the  measurement  time  tags  for  the  range  bias,  introducing  an  error  that  is 


3.3.  RESULTS 


19 


Table  3.1:  Relative  orbit  position  and  velocity  error  statistics  with  respect  to  laser  BERT  for  the 
Relative  State  Kalman  Filter  and  the  NASA/ JSC  filter  [4]. 


Position  (m) 

Velocity  (cm/s) 

Mean  STD  RMS 

Mean  STD  RMS 

RSKF 

Radial 

In-Track 

Cross-Track 

-2.8  5.8  6.4 
-4.2  3.2  5.3 
0.2  1.3  1.3 

-1.0  2.0  2.2 

-0.4  1.8  1.8 

0.1  1.3  1.3 

RSS 

8.4 

3.1 

JSC 

Radial 

In-Track 

Cross-Track 

-1.5  7.5  7.6 
-2.6  2.9  3.9 
-0.2  2.9  2.9 

-2.0  8.7  8.9 
0.0  6.1  6.1 
-1.2  4.0  4.2 

RSS 

9.0 

11.6 

most  likely  below  the  noise  level. 


3.3  Results 

Given  noisy,  asynchronous  observations  from  the  drag-laden,  maneuvering  Shuttle  and  the  relatively 
well-behaved  SPAS,  the  RSKF  obtains  excellent  results  considering  its  simplicity.  In  fact,  the  results 
are  comparable  or  better  than  results  obtained  in  a  NASA/JSC  study  detailed  in  [4].  Table  3.1 
shows  the  position  and  velocity  error  statistics  for  the  same  73  minute  period  using  the  same 
truth  reference  (the  laser  BERT)  for  the  RSKF  and  for  the  filter  developed  by  JSC.  The  JSC 
filter  estimates  the  relative  state  by  differencing  absolute  orbit  estimates  of  each  vehicle.  When  at 
least  four  common-view  satellites  are  visible,  the  JSC  filter  solves  for  the  orbits  with  only  those 
in  common  view,  using  the  same  filter-estimated  range  bias  for  a  given  GPS  SV .  If  less  than  four 
common  satellites  are  in  view,  all  channels  are  processed  for  both  receivers  [4], 

The  RSKF  results  are  obtained  with  the  following  initial  variances:  4000,  1000,  and  100  m2  for 
RIC  relative  position;  0.004,  0.004,  and  0.001  m2/s2  for  RIC  relative  velocity;  2000  m2  for  relative 
clock  bias;  and  25  m2/s2  for  relative  clock  frequency.  The  initial  spectral  densities  are  as  follows: 
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le-5,  4e-6,  8e-6  m2/s  for  RIC  relative  position;  3e-7,  4e-8,  8e-8  m2/s3  for  RIC  relative  velocity; 
0.01  m2/s  for  relative  clock  bias;  and  0.008  m2/s3  for  relative  clock  frequency.  The  noise  on  the 
measurement  is  estimated  to  be  approximately  3  m,  so  the  scalar  measurement  noise  value  R  is 
set  to  9  m2.  It  becomes  apparent  in  subsequent  plots  that  additional  tuning  must  be  performed  so 
that  the  filter  3<x  RMS  values  more  accurately  reflect  the  behavior  of  the  system. 

To  illustrate  the  nature  of  the  rendezvous,  Figs.  3.1,  3.2,  and  3.3  show  the  relative  position  and 
velocity  states.  The  in-plane  relative  trajectories  for  the  filter  and  laser  BERT  are  shown  in  Fig. 
3.1.  Note  that  motion  is  to  the  right  and  the  Earth  is  down.  It  can  be  seen  from  this  plot  that  the 
rendezvous  data  set  begins  at  a  range  of  approximately  1.8  km  and  ends  at  SPAS  capture.  Figs.  3.2 
and  3.3  show  all  the  relative  position  and  velocity  states  as  determined  by  the  RSKF  and  the  laser 
BERT.  These  plots  reveal  the  large  dynamic  motion  of  the  Shuttle  during  the  first  fifteen  minutes 
as  it  narrows  the  relative  range  from  1.8  km  to  200  m.  After  the  first  fifteen  minutes,  the  relative 
motion  becomes  more  benign,  but  is  still  difficult  to  model  because  the  Shuttle  is  maintaining  its 
position  along  the  R-bar.  That  is,  the  Shuttle  is  matching  the  velocity  of  the  SPAS  from  a  lower 
altitude. 

Another  indicator  of  the  large  motions  at  the  beginning  of  the  data  set  is  the  plot  of  the  delta-V 
magnitudes  in  Fig.  3.4.  For  the  purpose  of  velocity  aiding  in  the  RSKF,  a  range  of  minimum 
thresholds  was  evaluated.  A  threshold  of  1.5  cm/s  was  determined  to  provide  the  best  accuracy. 
The  figure  shows  a  line  at  1.5  cm/s,  with  the  points  above  it  marked  by  squares  to  indicate  those 
used  in  state  propagation. 

The  relative  position  errors  with  respect  to  the  laser  BERT  in  the  radial,  in-track,  and  cross¬ 
track  directions  are  shown  in  Fig.  3.5.  The  3cr  bounds  determined  by  the  RSKF  variances  are 
shown  as  the  dashed  lines.  Fig.  3.6  shows  the  relative  velocity  errors  with  respect  to  the  laser 
BERT.  Again,  the  3<r  bounds  determined  by  the  filter  variances  are  shown  as  the  dashed  lines. 
Compare  the  behavior  of  these  error  plots  with  the  motion  shown  in  Figs.  3.2  and  3.3.  Clearly, 
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the  RSKF  has  difficulty  obtaining  accurate  results  during  the  highly  dynamic  motion  of  the  first 
fifteen  minutes,  even  with  velocity  aiding.  After  that  period,  however,  the  filter  estimates  become 
less  erratic. 

Notice  the  biases  and  systematic  variations  in  the  relative  position  error  for  the  radial  and 
in-track  directions.  These  are  reflected  in  the  mean  and  RMS  position  error  statistics  in  Table 
3.1.  Possible  sources  of  error  include  inaccurate  GPS  antenna  CG  offsets,  time  tagging  errors,  or 
simply  unknown/unmodeled  biases.  The  fact  remains,  however,  that  these  biases  can  be  removed 
by  adjusting  the  body-fixed  antenna  CG  offsets  by  the  bias  amount.  Doing  so  lowers  the  position 
RMS  errors  to  the  same  level  as  the  standard  deviations,  thereby  reducing  the  total  RSS  position 
error  from  8.4  m  to  6.9  m.  There  is  only  a  minimal  improvement  in  the  velocity  statistics. 

To  ensure  that  the  RSKF  CW  propagation  is  not  introducing  a  bias  to  the  results,  a  hardware- 
in-the-loop  simulation  of  a  50  km  baseline  used  in  [1]  was  used  to  verify  its  accuracy.  Details  given 
in  the  next  chapter  show  unbiased  results.  The  results  do  show  a  an  error  due  to  relative  J2  for 
the  50  km  baseline,  but  that  should  not  be  an  issue  with  this  much  smaller  baseline.  Clearly,  the 
source  of  the  bias  in  the  Shuttle/SPAS  data  set  requires  further  investigation. 

The  most  dramatic  result  from  this  data  set  is  the  improvement  in  relative  velocity  estimates 
by  the  RSKF  versus  the  JSC  results.  The  RSKF  total  RSS  velocity  error  of  3.1  cm/s  is  a  73% 
improvement  over  the  11.6  cm/s  value  from  JSC  (see  Table  3.1).  Both  filters  include  velocity  aiding 
using  the  Shuttle  IMU  accelerometer  data.  In  fact,  both  filters  use  virtually  the  same  data,  so  the 
improvement  in  velocity  accuracy  is  most  likely  an  artifact  of  using  single  difference  observations 
in  the  RSKF. 

Finally,  the  plots  of  clock  behavior  are  given  in  Figs.  3.7  and  3.8.  Fig.  3.7  shows  the  estimates 
of  clock  bias  and  frequency  from  the  RSKF  and  from  GEODE,  the  absolute  orbit  determination 
utility  used  in  conjunction  with  the  RSKF.  The  difference  of  these  estimates  is  shown  in  Fig.  3.8. 
The  statistics  on  the  error  plot  reflect  values  after  the  first  five  minutes.  Also,  the  statistics  for  the 
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bias  error  include  extreme  outliers  that  occur  near  the  steps  of  1  ms  shown  in  Fig.  3.7.  If  outliers 
greater  than  100  /is  are  omitted  (21  of  1981  points),  the  RMS  error  drops  to  less  than  1  /is. 

These  two  plots  reflect  the  need  for  a  more  accurate  clock  model  in  the  RSKF.  Long  outages 
in  GPS  data  are  obvious  from  the  sections  of  constant  frequency  in  Fig.  3.7.  Corresponding  bias 
errors  can  be  seen  Fig.  3.8.  While  frequency  may  be  considered  constant  over  spans  of  seconds, 
the  assumption  quickly  becomes  invalid  for  outages  of  minutes.  More  appropriate  clock  models  are 
currently  being  evaluated. 


3.4  Conclusion 

In  a  relatively  simple  form,  the  RSKF  has  demonstrated  on  flight  data  the  ability  to  predict 
relative  position  and  velocity  to  accuracy  levels  comparable  to  or  better  than  the  absolute  state 
filter  developed  by  NAS  A/ JSC.  Further  improvements  in  positioning  accuracy  may  be  obtained 
by  identifying  the  source  of  the  relative  position  estimate  biases,  which  appear  in  both  the  RSKF 
and  the  JSC  filter  results  to  differing  degrees.  There  is  also  the  potential  for  additional  overall 
improvements  in  accuracy  by  implementing  a  relative  state  propagator  that  accounts  for  the  effects 
of  relative  gravity  and  drag  perturbations.  The  most  dramatic  improvement,  however,  is  gained 
by  using  more  precise  observations  and  techniques,  namely  carrier  phase  processing  coupled  with 
synchronization  of  the  vehicle  clocks.  The  potential  of  the  RSKF  with  synchronous  and  much  more 
accurate  measurements  is  shown  in  the  next  chapter. 


Vertical  (m) 
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Figure  3.2:  Radial,  in-track,  and  cross-track  relative  positions  for  the  filter  solutions  and  laser 
BERT. 
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Figure  3.3:  Radial,  in-track,  and  cross-track  relative  velocities  for  the  filter  solution  and  laser 
BERT. 
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Figure  3.4:  Shuttle  IMU-sensed  AV  magnitudes, 
signifying  the  points  used  in  state  propagation. 


Figure  3.6:  Relative  velocity  error  in  radial,  in-track,  and  cross-track  directions  with  respect  to 
laser.  Dotted  lines  show  filter  3 a  formal  error. 
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Figure  3.7:  Relative  clock  bias  and  frequency  estimates  from  GEODE  absolute  OD  utility  and  from 
RSKF.  Both  receivers  adjust  their  clocks  to  stay  within  0.5  ms  of  GPS  time.  Note  the  outages  in 
GPS  data  where  frequency  from  RSKF  becomes  constant. 
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Figure  3.8:  Relative  clock  bias  and  frequency  comparison  between  GEODE  and  RSKF  estimates. 
Statistics  do  not  include  first  5  minutes.  Bias  statistics  include  extreme  outliers  that  occur  near 
the  steps  (see  Fig.  3.7). 


Chapter  4 


Hardware  Simulation  Results 


Because  the  Shuttle/SPAS  data  include  large  maneuvers  and  noisy  data,  it  is  difficult  to  gauge  the 
potential  accuracy  of  the  RSKF  from  those  results.  Therefore,  I  ran  a  simplified  version  of  the 
RSKF  on  hardware  simulated  data  prepared  by  Patrick  Binning  at  NRL.  The  simulation  includes 
no  maneuvering,  so  the  maneuver  compensation  portion  of  the  RSKF  algorithm  is  not  utilized. 
Moreover,  Binning  provided  data  sets  pre-processed  through  the  NRL  orbit  determination  utility 
OCEAN  with  time  tags  already  aligned,  so  measurement  interpolation  is  not  necessary  in  the  RSKF . 
Finally,  no  antenna  CG  offsets  are  modeled,  so  that  feature  of  the  RSKF  is  also  not  utilized  here. 
What  remains  to  test  is  the  core  of  the  RSKF:  the  CW  propagation  model  and  the  measurement 
equations. 

The  details  of  the  hardware  simulations  can  be  found  in  [1].  The  results  shown  in  this  chapter 
are  from  only  one  of  the  many  simulations.  Binning  models  a  Front  and  Back  satellite  for  each 
simulation.  The  Back  satellite  follows  the  Front  with  along-track  separations  of  generally  5,  50, 
or  500  km.  The  scenario  given  here  is  for  a  separation  of  50  km,  with  initial  orbit  elements  given 
in  Table  4.1.  The  GPS  data  are  collected  from  an  Allen  Osborne  Associates  TurboRogue  receiver 
attached  to  a  Northern  Telecom  GPS  RF  signal  simulator  (an  earlier  version  of  the  simulator  to 


31 


32 


CHAPTER  4.  HARDWARE  SIMULATION  RESULTS 


Table  4.1:  Initial  orbital  elements  for  50-km  separation  hardware  simulation  performed  by  Binning 

ro¬ 


a 

7063.145  km 

e 

0.001 

i 

66.05  deg 

a 

0.0  deg 

U) 

0.0  deg 

M  (Front) 

270.162  deg 

M  (Back  50  km) 

269.756  deg 

Epoch  Date:  Sept.  20,  1996 
YUMA  Constellation:  GPS  Wk  871 

be  used  in  the  time  transfer  hardware  simulation).  Binning  set  every  controllable  error  source 
including  ionosphere  to  zero,  so  there  is  very  little  noise  on  the  measurements. 

The  observations  input  to  the  RSKF  are  LI  P-code  pseudoranges  output  from  OCEAN  with 
corrected  time  tags.  The  pseudoranges  for  both  receivers  are  also  interpolated  to  the  nearest  integer 
second  within  OCEAN,  allowing  explicit  differencing  in  the  RSKF.  Input  to  the  RSKF  also  includes 
the  OCEAN  ephemeris  for  the  target  vehicle  and  GPS  SV  positions  at  the  measurement  times.  Of 
the  9.5-hour  data  set,  I  processed  a  section  of  approximately  2.5  hours  that  has  good  common 
visibility  of  the  GPS  constellation. 


4.1  Results 

Table  4.2  shows  summary  statistics  of  the  results  as  compared  to  the  truth  output  from  the  sim¬ 
ulation.  The  statistics  do  not  include  the  first  20  minutes  during  which  the  filter  converges.  The 
corresponding  relative  position  and  velocity  error  plots  are  shown  in  Figs.  4.1  and  4.2.  Table  4.2 
shows  that  the  RSS  relative  position  error  is  only  26.1  cm  for  the  50  km  separation  distance.  On 
the  other  hand,  there  is  a  relatively  large  cross-track  velocity  error  that  causes  the  RSS  relative 
velocity  error  to  be  almost  5  cm/s.  While  the  in-plane  RMS  velocities  are  both  around  0.4  cm/s, 
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Table  4.2:  Relative  orbit  position  and  velocity  error  statistics  with  respect  to  simulation  truth  for 
the  Relative  State  Kalman  Filter  for  50  km  separation. 


Position  (cm) 

Velocity  (cm/s) 

Mean 

STD 

RMS 

Mean 

STD 

RMS 

Radial 

-4.5 

18.0 

18.6 

0.03 

0.41 

0.41 

In-Track 

-4.1 

12.3 

12.9 

0.02 

0.37 

0.37 

Cross-Track 

3.9 

12.4 

13.0 

-1.36 

4.56 

4.75 

RSS 

26.1 

4.78 

the  RSKF  appears  to  have  difficulty  resolving  the  relative  cross-track  motion. 

For  comparison  of  the  errors  to  their  relative  state  values,  Figs.  4.3  and  4.4  show  the  relative 
position  and  velocity,  respectively,  for  both  the  RSKF  (solid)  and  truth  (dotted).  Not  surprisingly, 
the  most  notable  feature  is  the  difference  in  cross-track  relative  velocity  estimates.  According  to  the 
RSKF,  there  is  essentially  no  cross-track  relative  motion,  while  the  truth  state  is  oscillating  with  the 
period  of  the  orbit  to  a  maximum  of  6  cm/s.  The  most  likely  cause  of  this  discrepancy  is  a  relative 
J2  effect.  For  a  50  km  along-track  separation  distance,  the  relative  J2  effect  will  be  the  difference 
of  small  accelerations  at  low  geocentric  latitudes,  and  the  difference  of  larger  accelerations  at  high 
geocentric  latitudes.  Correspondingly,  the  times  when  the  true  relative  cross-track  velocity  is  zero 
correspond  to  a  geocentric  latitude  of  zero;  whereas,  the  times  of  maximum  relative  cross-track 
velocity  correspond  to  maximum  geocentric  latitude. 

This  result  highlights  a  limitation  of  the  current  RSKF  implementation  with  a  simple  linear  CW 
model  and  no  velocity  dependence  in  the  measurement  equation  (Eq.  2.15).  This  effect  does  not 
show  up  in  the  Shuttle/SPAS  solution  because  of  the  close  proximity  of  the  vehicles  (within  500  m 
for  most  of  the  time).  Compensation  for  the  relative  J2  effect,  however,  should  improve  the  results 
on  data  sets  with  relatively  large  separations. 

Another  feature  to  note  in  Fig.  4.1  is  the  correlation  between  the  increase  in  formal  position 
error  and  the  decrease  in  common  GPS  satellites  in  view,  shown  in  Fig.  4.5.  The  number  of 
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common-view  SVs  drops  to  three  near  the  22  min,  45  min,  and  60  min  marks,  with  successively 
longer  durations.  Note  the  corresponding  increases  in  formal  position  error,  particularly  in  the 
radial  direction.  Other  portions  of  the  9.5-hour  data  set  include  worse  common  visibility  (including 
no  common  satellites),  at  which  point  the  solution  begins  to  diverge.  This  sensitivity  of  the  RSKF 
(or  OCEAN  absolute)  position  estimates  to  less  than  four  common-view  satellites  requires  further 
investigation. 

Fig.  4.6  shows  a  comparison  of  the  RSKF  relative  clock  estimates  to  the  difference  of  the 
OCEAN  absolute  estimates.  No  truth  reference  for  clock  bias  and  frequency  exists  for  this  hardware 
simulation,  so  the  comparison  is  made  to  results  from  another  filter.  The  relative  clock  estimates 
themselves  are  shown  in  Fig.  4.7.  Note  the  huge  relative  frequency  step  in  Fig.  4.7.  Due  to  the 
existence  of  that  step,  the  process  noise  for  that  state  must  be  large,  as  indicated  by  the  large  formal 
error  in  the  frequency  error  plot  of  Fig.  4.6.  Also,  the  factor  of  40  step  causes  range  residuals  so 
large  that  the  RSKF  residual  editor  never  recovers.  Therefore,  the  residual  editing  feature  is  not 
implemented  for  this  data  set.  This  unexpected  frequency  behavior  highlights  another  area  for 
improvement  in  the  RSKF.  The  residual  editing  scheme  must  be  enhanced  to  allow  residual  editing 
while  suffering  only  short  outages  due  to  anomalous  behavior. 


4.2  Conclusion 

The  low  noise  and  large  separation  of  this  simulated  data  allows  better  evaluation  of  the  filter 
algorithms  than  the  Shuttle/SPAS  flight  data.  While  the  results  indicate  a  couple  shortcomings  of 
the  algorithms,  what  is  also  clear  is  that  the  RSKF  in  its  current  state  can  resolve  relative  positions 
to  less  than  20  cm  in  each  direction  over  a  50  km  baseline.  Furthermore,  with  compensation  for 
relative  J2,  these  results  indicate  that  relative  velocity  could  be  resolved  to  0.4  cm/s  or  better 
in  each  direction.  Carrier  phase  processing  should  improve  the  results  even  further,  though  not 
dramatically  because  of  the  low  noise  on  the  P-code  pseudoranges.  In  summary,  the  RSKF  provides 
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accurate  results  for  most  states,  but  requires  modification  to  improve  the  accuracy  of  the  cross-track 
relative  velocity. 
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Figure  4.1:  Radial,  in-track,  and  cross- track  relative  position  errors  for  the  RSKF  with  respect  to 
simulation  truth.  Dotted  lines  show  filter  la  formal  error. 
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Figure  4.2:  Radial,  in-track,  and  cross-track  relative  velocity  errors  for  the  RSKF  with  respect  to 
simulation  truth.  Dotted  lines  show  filter  lcr  formal  error. 
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Figure  4.3:  Relative  position  in  radial,  in-track,  and  cross-track  directions  for  RSKF  and  truth. 
Note  that  the  scale  on  the  in-track  ordinate  is  km. 
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Figure  4.4:  Relative  velocity  in  radial,  in-track,  and  cross-track  directions  for  RSKF  and  truth. 
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Figure  4.6:  Relative  clock  bias  and  frequency  errors  for  RSKF  with  respect  to  difference  of  absolute 
clock  states  from  NRL  OCEAN  estimates.  Dotted  lines  show  filter  lcr  formal  error.  Mean  clock 
bias  error  is  subtracted  out. 


Chapter  5 


Time  Transfer  Hardware  Simulation 


The  algorithms  for  on-orbit  relative  navigation,  timing,  and  frequency  will  be  tested  primarily  using 
newly-generated  hardware-in-the-loop  simulations  which  will  include  a  truth  reference  for  the  clock 
states.  The  following  sections  describe  the  hardware  to  be  used  in  the  test  bed  and  the  procedures 
for  testing. 


5.1  Test  Bed 

Fig.  5.1  represents  the  basic  test  bed  for  two  GPS  receivers  being  tested  in  parallel.  Two  Global 
Simulation  Systems  (GSS)  GPS  signal  simulators  will  provide  simultaneous  signals  to  separate 
Trimble  Advanced  Navigation  Sensor  (TANS)  six-channel  GPS  receivers.  The  two  receivers  will 
be  modeled  in  neighboring  orbits.  The  1  pulse-per-second  (PPS)  output  from  each  receiver  will 
be  measured  in  both  an  absolute  and  a  relative  sense  by  interval  counters.  A  hydrogen  maser 
operating  at  5  MHz  will  provide  an  accurate  reference  oscillator  signal  to  the  simulators  as  well 
as  to  the  interval  counters.  The  maser  will  also  provide  a  reference  1  PPS  signal  to  the  interval 
counter  measuring  the  absolute  accuracy  of  the  receiver  clocks.  Knowledge  of  the  absolute  clocks 


43 


44 


CHAPTER  5.  TIME  TRANSFER  HARDWARE  SIMULATION 


Data  for 
Sat  2 


Figure  5.1:  Time  transfer  test  bed.  Absolute  and  relative  measurement  of  the  1  PPS  out  of  the 
GPS  receivers  is  actually  performed  in  separate  runs  of  identical  simulations. 


is  desired  to  help  identify  any  anomalous  relative  behavior. 

The  GSS  simulators  produce  a  comprehensive  development  and  testing  environment.  The  simu¬ 
lator  consists  of  Multichannel  Satellite  Signal  Source  (MSSS)  and  a  host  computer.  The  MSSS  is  the 
actual  hardware  unit  that  generates  the  RF  signals  sent  to  the  receiver.  The  signal  generator  can  be 
considered  a  pseudorange-to-RF  converter  [20].  The  host  computer  drives  the  MSSS  and  processes 
the  results  of  the  simulation.  The  GSS  software  allows  specification  or  adjustment  of  almost  (if  not) 
all  parameters  of  the  simulation.  Some  examples  are  spacecraft  characteristics,  dynamics  limits, 
antenna  details,  orbit  initial  conditions  or  ephemeris,  maneuvers,  attitudes,  ionosphere,  multipath, 
and  GPS  constellation  configurations  and  signal  strengths.  NRL  is  currently  running  version  6.46 
of  the  GSS  firmware,  versus  version  1.20  when  Binning  performed  his  simulations  [1]. 


The  Trimble  receiver  is  a  six-channel,  Standard  Positioning  Service  (SPS)  receiver  which  uses 
the  C/A  code  on  LI.  The  stated  PPS  output  has  a  width  of  1  microsecond  with  the  falling  edge 
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synchronized  to  UTC  within  1  microsecond  [21].  In  a  conversation  with  Joe  White  of  NRL,  he 
indicated  that  the  resolution  of  the  TANS  PPS  should  be  better,  approximately  100  nanoseconds 
absolute.  Although  not  as  precise  as  desired,  the  benefit  of  using  the  TANS  is  the  fact  that  NRL 
has  two.  This  allows  synchronous  simulation  as  portrayed  in  Fig.  5.1.  If  only  one  receiver  is 
used,  separate  tests  must  be  performed,  making  synchronization  and  measurement  of  relative  clock 
behavior  problematic.  Also,  if  two  different  receiver  types  are  used,  that  requires  processing  two 
different  receiver  output  formats. 

Future  tests,  however,  may  include  the  Motorola  GPSR  dual-frequency,  twelve-channel  receiver 
developed  for  NRL.  These  tests  may  be  performed  in  series  or  in  parallel  with  a  TANS,  at  which 
point  the  limitations  mentioned  above  must  be  overcome.  The  outcome  of  preliminary  tests  with 
the  TANS  receivers  will  help  determine  the  direction  of  future  testing. 


5.2  Test  Plan 

The  test  bed  depicted  in  Fig.  5.1  actually  comprises  two  simulation  runs.  Due  to  potential 
impedance  problems  caused  by  splitting  the  PPS  output  from  the  GPS  receivers,  the  exact  same 
simulation  scenario  will  be  run  twice,  once  to  measure  the  absolute  clock  behavior  and  again  to 
measure  the  relative  clock  behavior.  The  relative  clock  behavior  cannot  be  determined  accurately 
from  the  absolute  measurements  because  they  are  not  taken  simultaneously.  The  interval  counter 
samples  the  channels  of  the  switch  unit  only  every  5  minutes.  However,  with  the  PPS  output  from 
the  receivers  plugged  directly  into  the  interval  counter  during  the  second  run,  measurements  of 
relative  clock  behavior  can  be  taken  at  1  Hz. 

Simulations  with  varying  baselines  will  be  run.  Because  the  TANS  receiver  can  see  only  a 
maximum  of  six  satellites,  separation  distances  greater  than  several  hundred  kilometers  may  not 
include  a  sufficient  number  of  common-view  GPS  satellites.  The  first  tests  will  include  separation 
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distances  of  10  km  and  100  km  with  both  along-track  and  cross-track  initial  offsets.  All  specifiable 
error  sources  will  be  set  to  zero  for  the  initial  tests  so  that  the  C/A  code  pseudoranges  from  the 
TANS  are  as  noiseless  as  possible.  The  measured  relative  clock  bias  and  frequency  will  be  evaluated 
for  sufficient  precision  before  proceeding.  Assuming  small-enough  time  resolution,  the  data  will  be 
processed  with  the  RSKF,  and  the  filter  models  and  algorithms  will  be  adjusted  to  achieve  the 
highest  possible  precision. 

If  the  TANS  proves  unreliable,  the  most  likely  course  of  action  is  to  employ  the  Motorola  GPSR 
receiver  in  series  simulations.  The  task  of  obtaining  an  accurate  truth  reference  for  the  relative 
clock  behavior  will  be  worked  out  with  the  help  of  committee  members  and  the  timing  experts  at 
NRL.  If  neither  GPS  space  receiver  can  resolve  time  to  the  desired  precision,  then  my  research  will 
attempt  to  estimate  the  relative  clock  states  to  the  measured  precision. 


5.3  Conclusions 

Successful  completion  of  this  research  will  produce  the  first  demonstration  of  precise  relative  nav¬ 
igation  and  timing  of  satellites  in  orbit  in  real  time  using  actual  hardware  data.  This  capability 
will  permit  applications  such  as  synchronization  of  observations  from  a  distributed  system  of  sparse 
aperture  radar  satellites.  Another  contribution  will  be  the  design  of  the  test  bed  for  simultaneous 
hardware  simulation  of  co-orbiting  vehicles  and  measurement  of  the  absolute  and  relative  clock 
behavior.  This  will  be  a  useful  tool  for  analyzing  the  relative  navigation  and  timing  performance  of 
GPS  receivers  for  future  formation  flying  architectures.  Overall,  this  will  be  an  interesting  research 
topic  with  the  potential  to  produce  valuable  new  tools  for  satellite-to-satellite  navigation  and  time 


transfer. 
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